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Abstract 

DDSCAT 7.1 is a freely available open-source Fortran-90 software package applying the "discrete 
dipole approximation" (DDA) to calculate scattering and absorption of electromagnetic waves by targets 
with arbitrary geometries and complex refractive index. The targets may be isolated entities (e.g., dust 
particles), but may also be 1-d or 2-d periodic arrays of "target unit cells", which can be used to study 
absorption, scattering, and electric fields around arrays of nanostructures. 

The DDA approximates the target by an array of polarizable points. The theory of the DDA and its 
implementation in DDSCAT is presented in Draine (1988) and Draine & Flatau (1994), and its extension 
to periodic structures (and near-field calculations) in Draine & Flatau (2008). DDSCAT 7.1 allows 
accurate calculations of electromagnetic scattering from targets with "size parameters" 27ra e ff / A ^ 25 
provided the refractive index m is not large compared to unity (|m — 1| ^ 2). DDSCAT 7.1 includes 
support for MP I, OpenMP, and the Intel® Math Kernel Library (MKL). 

DDSCAT supports calculations for a variety of target geometries (e.g., ellipsoids, regular tetrahedra, 
rectangular solids, finite cylinders, hexagonal prisms, etc.). Target materials may be both inhomoge- 
neous and anisotropic. It is straightforward for the user to "import" arbitrary target geometries into the 
code. DDSCAT automatically calculates total cross sections for absorption and scattering and selected 
elements of the Mueller scattering intensity matrix for specified orientation of the target relative to the 
incident wave, and for specified scattering directions. DDSCAT 7.1 can calculate scattering and absorp- 
tion by targets that are periodic in one or two dimensions. DDSCAT 7.1 also supports postprocessing 
to calculate E and B at user-selected locations in or near the target; a Fortran-90 code DDfleld for this 
purpose is included in the disnibution. 

This User Guide explains how to use DDSCAT 7.1 (release 7.1.0) to carry out electromagnetic 
scattering calculations. 
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1 Introduction 

DDSCAT is a software package to calculate scattering and absorption of electromagnetic waves by tar- 
gets with arbitrary geometries using the "discrete dipole approximation" (DDA). In this approximation 
the target is replaced by an array of point dipoles (or, more precisely, polarizable points); the electromag- 
netic scattering problem for an incident periodic wave interacting with this array of point dipoles is then 
solved essentially exactly. The DDA (sometimes referred to as the "coupled dipole approximation") was 
apparently first proposed by Purcell & Pennypacker (1973). DDA theory was reviewed and developed 
further by Draine (1988), Draine & Goodman (1993), reviewed by Draine & Flatau (1994) and Draine 
(2000), and recently extended to periodic structures by Draine & Flatau (2008). 

DDSCAT 7.1, the current release of DDSCAT, is an open-source Fortran 90 implementation of the 
DDA developed by the authorsQ It extends DDSCAT to treat structures that are periodic in one or two 
dimensions, using methods described by Draine & Flatau (2008). 

DDSCAT is intended to be a versatile tool, suitable for a wide variety of applications including 
studies of interstellar dust, atmospheric aerosols, blood cells, marine microorganisms, and nanostructure 
arrays. As provided, DDSCAT 7.1 should be usable for many applications without modification, but the 
program is written in a modular form, so that modifications, if required, should be fairly straightforward. 

The authors make this code openly available to others, in the hope that it will prove a useful tool. We 
ask only that: 

• If you publish results obtained using DDSCAT, please acknowledge the source of the code, and 
cite relevant papers, such as Draine & Flatau (1994), Draine & Flatau (2008), and this UserGuide 
(Draine & Flatau 2010). 

• If you discover any errors in the code or documentation, please promptly communicate them to the 
authors. 

• You comply with the "copyleft" agreement (more formally, the GNU General Public License) of 
the Free Software Foundation: you may copy, distribute, and/or modify the software identified as 
coming under this agreement. If you distribute copies of this software, you must give the recipients 
all the rights which you have. See the file doc/copylef t . txt distributed with the DDSCAT 
software. 

We also strongly encourage you to send email to draine@astro.princeton.edu identifying 
yourself as a user of DDSCAT; this will enable the authors to notify you of any bugs, corrections, or 
improvements in DDSCAT. Up-to-date information on DDSCAT can be found at 

|http : //www .astro . princeton . edu/~draine/DDSCAT . html| 
Information is also available at |http : //ddscat . wikidot . com/start| 

The current version, DDSCAT 7.1, uses the DDA formulae from Draine (1988), with dipole polariz- 
abilities determined from the Lattice Dispersion Relation (Draine & Goodman 1993, Gutkowicz-Krusin 
& Draine 2004). The code incorporates Fast Fourier Transform (FFT) methods (Goodman, Draine, & 
Flatau 1991). DDSCAT 7.1 includes capability to calculate scattering and absorption by targets that are 

1 The release history of DDSCAT is as follows: 

• DDSCAT 4b: Released 1993 March 12 

• DDSCAT 4bl: Released 1993 July 9 

• DDSCAT 4c: Although never announced, DDSCAT . 4 c was made available to a number of interested users beginning 1994 
December 18 

• DDSCAT 5a7: Released 1996 

• DDSCAT 5a8: Released 1997 April 24 

• DDSCAT 5a9: Released 1998 December 15 

• DDSCAT 5al0: Released 2000 June 15 

• DDSCAT 6.0: Released 2003 September 2 

• DDSCAT 6.1: Released 2004 September 10 

• DDSCAT 7.0: Released 2008 September 1 

• DDSCAT 7.1: Released 2010 February 7 
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periodic in one or two dimensions - arrays of nanostructures, for example. The theoretical basis for 
application of the DDA to periodic structures is developed in Draine & Flatau (2008). 

We refer you to the list of references at the end of this document for discussions of the theory and 
accuracy of the DDA [in particular, reviews by Draine and Flatau (1994) and Draine (2000), and recent 
extension to 1-d and 2-d arrays by Draine & Flatau (2008).]. In Sj4]we describe the principal changes 
between DDSCAT 7.1 and the previous releases. The succeeding sections contain instructions for: 

• compiling and linking the code; 

• running a sample calculation; 

• understanding the output from the sample calculation; 

• modifying the parameter file to do your desired calculations; 

• specifying target orientation; 

• changing the DIMENSIONing of the source code to accommodate your desired calculations. 

The instructions for compiling, linking, and running will be appropriate for a Linux system; slight 
changes will be necessary for non-Linux sites, but they are quite minor and should present no difficulty. 

Finally, the current version of this User Guide can be obtained from 
|http ://arxiv.org/abs/1002.1505| - you will be offered the options of downloading either 
Postscript or PDF versions of the User Guide. 



2 Applicability of the DDA 

The principal advantage of the DDA is that it is completely flexible regarding the geometry of the target, 
being limited only by the need to use an interdipole separation d small compared to (1) any structural 
lengths in the target, and (2) the wavelength A. Numerical studies (Draine & Goodman 1993; Draine & 
Flatau 1994; Draine 2000) indicate that the second criterion is adequately satisfied if 

\m\kd<\ , (1) 

where m is the complex refractive index of the target material, and k = 2ir/\, where A is the wavelength 
invacuo. This criterion is valid provide |m— 1| ;$ 3orso. When Im(m) becomes large, the DDA solution 
tends to overestimate the absorption cross section C a b s , and it may be necessary to use interdipole 
separations d smaller than indicated by eq. (H) to reduce the errors in C a b s to acceptable values. 

If accurate calculations of the scattering phase function (e.g., radar or lidar cross sections) are desired, 
a more conservative criterion 

|m|fcd<0.5 (2) 

will usually ensure that differential scattering cross sections dC sca /dQ, are accurate to within a few 
percent of the average differential scattering cross section C S ca/47r (see Draine 2000). 

Let V be the actual volume of solid material in the target|j If the target is represented by an array of 
N dipoles, located on a cubic lattice with lattice spacing d, then 

U = A^ 3 . (3) 

We characterize the size of the target by the "effective radius" 

a cff = (31//47T) 1 / 3 , (4) 

the radius of an equal volume sphere. A given scattering problem is then characterized by the dimen- 
sionless "size parameter" 

- , 27ra ff 
x = fca cff = — - — . (5) 



2 In the case of an infinite periodic target, V is the volue of solid material in one "Target Unit Cell". 
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The size parameter can be related to N and |m|fc<i: 



ss 2^ = 62,04 (NV'* _ _ (6) 



A |m| V !0 6 . 

Equivalently, the target size can be written 

A /A^ 1/3 



a ° ff = 9 - 873 H Uv ■ |m|fcd • (?) 

Practical considerations of CPU speed and computer memory currently available on scientific worksta- 
tions typically limit the number of dipoles employed to TV < 10 6 (see |T6]for limitations on TV due to 
available RAM); for a given N, the limitations on |m|fcd translate into limitations on the ratio of target 
size to wavelength. 

For calculations of total cross sections C a bs and C sca , we require |m|fcd < 1: 

A / N A 1/3 62.04 / N N 1/3 



a eff < 9.88— — - or x < — — — . (8) 

\m\ \10 6 J \m\ \10 6 J 

For scattering phase function calculations, we require \m\kd < 0.5: 

tnl A ( N A 1/3 31.02 / N \ 1/3 

a e ff < 4.94- — r — - or x < - — - — - . (9) 

\m\ VlO 6 J \m\ VlO 6 J 

It is therefore clear that the DDA is not suitable for very large values of the size parameter x, or very 
large values of the refractive index m. The primary utility of the DDA is for scattering by dielectric 
targets with sizes comparable to the wavelength. As discussed by Draine & Goodman (1993), Draine 
& Flatau (1994), and Draine (2000), total cross sections calculated with the DDA are accurate to a few 
percent provided N > 10 4 dipoles are used, criterion (Q]) is satisfied, and the refractive index is not too 
large. 

For fixed |m|fcd, the accuracy of the approximation degrades with increasing \m — 1|, for reasons 
having to do with the surface polarization of the target as discussed by Collinge & Draine (2004). With 
the present code, good accuracy can be achieved for \m — 1| < 2. 

Examples illustrating the accuracy of the DDA are shown in Figs.[T]-[2]which show overall scattering 
and absorption efficiencies as a function of wavelength for different discrete dipole approximations to a 
sphere, with N ranging from 304 to 59728. The DDA calculations assumed radiation incident along the 
(1,1,1) direction in the "target frame". Figs.[3]43show the scattering properties calculated with the DDA 
for x = ka = 7. Additional examples can be found in Draine & Flatau (1994) and Draine (2000). 

As discussed below, DDSCAT 7.1 can also calculate scattering and absorption by targets that are 
periodic in one or two directions - for examples, see Draine & Flatau (2008). 



3 DDSCAT 7.1 

3.1 What Does It Calculate? 

3.1.1 Absorption and Scattering by Finite Targets 

DDSCAT 7.1, like previous versions of DDSCAT, solves the problem of scattering and absorption by 
a finite target, represented by an array of polarizable point dipoles, interacting with a monochromatic 
plane wave incident from infinity. DDSCAT 7.1 has the capability of automatically generating dipole 
array representations for a variety of target geometries (see §1% and can also accept dipole array repre- 
sentations of targets supplied by the user (although the dipoles must be located on a cubic lattice). The 
incident plane wave can have arbitrary elliptical polarization (see £ j22l) . and the target can be arbitrarily 
oriented relative to the incident radiation (see CPTI i. The following quantities are calculated by DDSCAT 
7.1 : 
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1 2 3 4 5 6 7 8 9 10 11 18 13 
x=ka 



Figure 1: Scattering and absorption for a sphere with m = 1.33 + O.Oli. The upper panel shows the exact values of Q sca and 
Qabs> obtained with Mie theory, as functions of x — ka. The middle and lower panels show fractional errors in Q SC a and Q a bs, 
obtained using DDSCAT with polarizabilities obtained from the Lattice Dispersion Relation, and labelled by the number N of 
dipoles in each pseudosphere. After Fig. 1 of Draine & Flatau (1994). 
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x = ka 



Figure 2: Same as Fig.fU but for m = 2 + i. After Fig. 2 of Draine & Flatau (1994). 



3 DDSCAT7.1 



10 




Figure 3: Differential scattering cross section for m — 1.33 + O.Oli pseudosphere and ka = 7. Lower panel shows fractional 
error compared to exact Mie theory result. The N = 17904 pseudosphere has \m\kd — 0.57, and an rms fractional error in da/dQ 
of 2.4%. After Fig. 5 of Draine & Flatau (1994). 




Figure 4: Same as Fig.[5]but for m = 2 + i. The N = 59728 pseudosphere has \m\kd = 0.65, and an rms fractional error in 
da/dQ of 6.7%. After Fig. 8 of Draine & Flatau (1994). 
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• Absorption efficiency factor Q a bs = C a b s /7nXg ff , wnere Cabs is the absorption cross section; 

• Scattering efficiency factor <5 sca = C sca /ira^ s , where C sca is the scattering cross section; 

• Extinction efficiency factor Q CKt = Q sca + Qab s ; 

• Phase lag efficiency factor Q p h a , defined so that the phase-lag (in radians) of a plane wave after 
propagating a distance L is just n t Q p ha,^a^sL, where n t is the number density of targets. 

• The 4x4 Mueller scattering intensity matrix Sij describing the complete scattering properties of 
the target for scattering directions specified by the user (see ^24b . 

• Radiation force efficiency vector Q la( j (see 915) . 

• Radiation torque efficiency vector Qr (see Ci"5l ). 

In addition, the user can choose to have DDSCAT 7.1 store the solution for post-processing. A separate 
Fortran-90 program DDfleld (see $T7\ is provided that can readily calculate E and B at user-selected 
locations. 



3.1.2 Absorption and Scattering by Periodic Arrays of Finite Structures 

DDSCAT 7.1 includes the capability to solve the problem of scattering and absorption by an infinite 
target consisting of a 1-d or 2-d periodic array of finite structures, illuminated by an incident plane wave. 
The finite structures are themselves represented by arrays of point dipoles. 

The electromagnetic scattering problem is formulated by Draine & Flatau (2008), who show how the 
problem can be reduced to a finite system of linear equations, and solved by the same methods used for 
scattering by finite targets. 

The far- field scattering properties of the 1-d and 2-d periodic arrays can be conveniently represented 
by a generalization of the Mueller scattering matrix to the 1-d or 2-d periodic geometry - see Draine & 
Flatau (2008) for definition of S^ d) (M, () and S^ d) (M, N). For targets with 1-d periodicity, DDSCAT 
7.1 calculates S^ d) {M, () for user-specified M and £. For targets with 2-d periodicity, DDSCAT 7.1 

calculates s!-j (M, N) for both transmission and reflection, for user-specified (M, N). 

As for finite targets, the user can choose to have DDSCAT 7.1 store the calculated polarization 
field for post-processing. A separate Fortran-90 program DDfleld is provided to calculate E and B at 
user-selected locations, which can be located anywhere (including within the target). 



3.2 Application to Targets in Dielectric Media 

Let oj be the angular frequency of the incident radiation. For many applications, the target is essen- 
tially in vacuo, in which case the dielectric function e which the user should supply to DDSCAT is the 
actual complex dielectric function etarget(w), or complex refractive index retarget (oj) = Retarget of the 
target material. 

However, for many applications of interest (e.g., marine optics, or biological optics) the "target" 
body is embedded in a (nonabsorbing) dielectric medium (e.g., water at visible wavelengths) with (real) 
dielectric function e met a um (uj), or (real) refractive index m me dium(w) = ^medium- DDSCAT is fully 
applicable to these scattering problems, except that: 

• The complex "dielectric function" or complex "refractive index" supplied to DDSCAT should be 
the relative dielectric function 

e ( w ) = target H (10) 

^medium ) 

or relative refractive index: 

m{u)= m target H (U) 

m-modium(w) 

• The wavelength A specified in ddscat . par should be the wavelength in the medium: 

A = Avac , (12) 



'■medium 



where A vac = 2ttc/uj is the wavelength in vacuo. 
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Because DDSCAT will be using the relative refractive index m(u>) and the wavelength in the medium, 
the file supplying the refractive index as a function of wavelength should give the relative refractive 
index or relative dielectric function as a function of wavelength in the medium. Those using DDSCAT 
for targets that are not in vacuo must be attentive to this when preparing the input table(s) giving the 
relative dielectric function or refractive index of the target material(s). 

The absorption, scattering, extinction, and phase lag efficiency factors Q a b s , Qsca, and Q cx t cal- 
culated by DDSCAT will then be equal to the physical cross sections for absorption, scattering, and 
extinction divided by 7ro^ ff . For example, the attenuation coefficient for radiation propagating through a 
medium with a number density n t of scatterers will be just a = n t Q cy r t Tra^. s . Similarly, the phase lag 
(in radians) after propagating a distance L will be ntQp\ l& 'Ka'^ s L. 

The elements Sij of the 4x4 Mueller scattering matrix S calculated by DDSCAT for finite targets 
will be correct for scattering in the medium: 



= (—)' 

V27rr J 



S • Ii„, (13) 



where Ij n and I sca are the Stokes vectors for the incident and scattered light (in the medium), r is the 
distance from the target, and A is the wavelength in the medium (eq.[T2l). See S|24]for a detailed discussion 
of the Mueller scattering matrix. 

The time-averaged radiative force and torque (see §15[ on a finite target in a dielectric medium are 

F ra d = QprTTdcff^rad , (14) 
Trad = QrTTacffWrad^- , (15) 

where the time-averaged energy density is 

El 

^rad — ^medium ~Z ; (16) 

where Eq cos(ujt + <fi) is the electric field of the incident plane wave in the medium. 



4 What's New? 

DDSCAT 7.1 differs from DDSCAT 6.1 in a number of ways, including: 

1. N.B.: The structure of the parameter file ddscat . par has been changed: ddscat . par files 
that were used with DDSCAT 6.1 will need to be modified. See Sj8]and AppendixlAl 

2. DDSCAT 7.0 can calculate scattering by targets consisting of anisotropic materials with arbitrary 
orientation. Two new target options ANIFRMFIL and SPH_ANI_N are provided that make use of 
this capability. 

ANIFRMFIL reads an input file with a list of dipoles for a general anisotropic target. For each 
dipole, the input file gives the dipole location, composition for each of three principal axes of 
the local dielectric tensor, and the orientation of the local "Dielectric Frame" (DF) relative to the 
"Target Frame" (TF). 

Target option SPH_ANI_N creates target consisting of N anisotropic spheres. An input file lists, 
for each sphere, the location, radius, composition corresponding to each of three principal axes of 
dielectric tensor, and the orientation of the DF for this sphere relative to the TF. 

3. The user now has the option of specifying the scattering directions either relative to the Lab Frame 
(as in DDSCAT 6. 1) or relative to the Target Frame. This is selected by providing a new character 
variable, either L FRAME or TFRAME in the parameter file ddscat .par. Note that this choice 
must be specified in ddscat . par. 

4. The user must now specify the value of parameter IWRPOL in the parameter file ddscat . par. 
See Sj8]and Appendix|A] 
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5. A separate Fortran-90 code DDfleld is provided that can read the stored polarization array that is 
written if IWRPOL = 1 and use this data to calculate E and B at user-specified positions in or near 
the target. 

6. DDSCAT 7.0 can solve the problem of scattering and absorption by a target that is periodic in 
one dimension (and finite in the other two dimensions) or periodic in two dimensions (and finite 
in the third) and illuminated by a monochromatic plane wave, using periodic boundary condi- 
tions (PBC). The theory is described in Draine & Flatau (2008). A number of built-in target op- 
tions are provided (see S|20]l, including BISLINPBC, CYLNDRPBC, DSKLYRPBC, HEXGONPBC, 
LYRSLBPBC, RCTGL_PBC, SPHRN_PBC, 

7. DDSCAT 7.0 uses dynamic memory allocation. An additional line is now required in the dds cat .par 
file to inform DDSCAT 7.1 how much memory should initially be allocated to safely generate the 
requested target. Once the target has been generated, DDSCAT 7.1 reevaluates the memory re- 
quirements, and then allocates only the memory required to carry out the scattering calculation for 

the chosen target. 

8. The user can choose to compile DDSCAT 7.0 in either single- or double-precision. 

9. DDSCAT 7.0 includes built-in support for OpenMP; if this is enabled, certain parts of the calcula- 
tion can be distributed over available cpu cores. This can speed execution on multi-core systems. 

10. DDSCAT 7.0 includes built-in support to use the Intel® Math Kernel Library (MKL) for calcula- 
tion of FFTs. 

1 1 . DDSCAT 7.1 is distributed with a set of examples, including both dds cat . par files and selected 
output files. These should be useful to new users. 

12. DDSCAT 7.1 (version 7.1.0) allows the user to study periodic targets with the target unit cell 
(TUC) geometry (i.e., dipole locations) and composition provided via an input file [new target 
options FRMFILPBC ( ^20. lb and ANIFILPBC ( ^20.2H . 

13. DDSCAT 7.1 (version 7.1.0) supports new target options SLAB_HOLE (rectangular block with 
cylindrical hole) and SLBHOLPBC periodic array of rectangular blocks with cylindrical holes. 
SLBHOLPBC can be used to create a uniform slab with uniformly-spaced cylindrical holes in a 
rectangular pattern. 

14. DDSCAT 7.1 (version 7.1.0) includes a new version of DDfleld, with modification to the structure 
of the input file ddf ie Id . in to facilitate sampling the E and B fields at many points. 

5 Obtaining and Installing the Source Code 

DDSCAT 7.1 is written in standard Fortran-90 with a single extension: it uses the Fortran-95 standard 
library call CPU.TIME to obtain timing information. DDSCAT 7.1 is therefore portable to any system 
having a f90 or f95 compiler. It has been successfully compiled with many different compilers, including 
gfortran, g95, ifort, pgf 77, andNAG®f 95. 

It is possible to use DDSCAT 7.1 on PCs running Microsoft Windows operating systems, including 
Vista and Windows 7. More information on how to do this can be found at 
|http : //code . google . com/p/ddscat] 

The remainder of this section will assume that the installation is taking place on a Unix, Linux or Mac 
OSX system with the standard developer tools (e.g., tar, make, and a f90 or f95 compiler) installed. 
The complete source code for DDSCAT 7.1 is provided in a single gzipped tarfile. 
There are 4 ways to obtain the source code: 

• Point your browser to the DDSCAT home page: 

[http : / /www .astro . princeton . edu /^draine /DP SCAT . html| 
and right=click on the link to ddscat7 .1.0. tgz . 

• Point your browser to 

|http : / /code . google . com/p/ddscat71 and download the latest release. 
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• Point your browser to 

[ftp : / / ftp . astro .princeton . edu/ draine/scat/ddscat/7 . 1 . 0| 
right-click on ddscat7 .1.0 . tgz . 

• Use anonymous ftp to ftp.astro.princeton.edu. 
In response to "Name:" enter anonymous. 

In response to "Password:" enter your full email address. 

cd draine/scat/ddscat/ver7 . 1 . 
binary ( to enable transfer of binary data) 
get ddscat7 . 1 . . tgz 

After downloading ddscat7 .1.0. tgz into the directory where you would like DDSCAT 7.1 to re- 
side (you should have at least 10 Mbytes of disk space available), the source code can be installed as 
follows: 

If you are on a Linux system, you should be able to type 

tar xvzf ddscat7 . 1 . . tgz 
which will"extract" the files from the gzipped tarfile. If your version of "tar" doesn't support the "z" 
option (e.g., you are running Solaris) then try 

zcat ddscat7 . 1 . . tgz | tar xvf - 
If neither of the above work on your system, try the two-stage procedure 

gunzip ddscat7 . 1 . . tgz 

tar xvf ddscat7 . 1 . . tar 
A disadvantage of this two-stage procedure is that it uses more disk space, since after the second step 
you will have the uncompressed tarfile ddscat 7 . 1 . . tar - about 3.8 Mbytes - in addition to all the 
files you have extracted from the tarfile - another 4.6 Mbytes. 

Any of the above approaches should create subdirectories src, doc, and examples.exp. The 
source code will be in subdirectory src, and documentation in subdirectory doc. The examples.exp 
subdirectory contains sample ddscat .par files as well as output files for various target options, in- 
cluding both isolated targets and infinite periodic targets. 

6 Compiling and Linking 

In the discussion below, it is assumed that the DDSCAT 7.1 source code has been installed in a directory 
DDA/ src. The instructions below assume that you are on a Unix, Linux, or Mac OSX system. 

6.1 The default Makefile 

It is assumed that the system has the following already installed: 

• a Fortran-90 compiler (e.g., gf ortran, g95, Intel® if ort, or NAG® f 95) . 

• cpp - the "C preprocessor". 

There are a number of different ways to create an executable, depending on what options the user wants: 

• what compiler and compiler flags? 

• single- or double-precision? 

• enable OpenMP? 

• enable MKL? 

• enable MPI? 

Each of the above choices needs requires setting of appropriate "flags" in the Makefile. 
The default Makefile has the following "vanilla" settings: 

• gfortran -02 

• single-precision arithmetic 

• OpenMP not used 
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• MKL not used 

• MPI not used. 

To compile the code with the default settings, simply position yourself in the directory DDA/src, and 
type 

make ddscat 

If you have gf ortran and cpp installed on your system, the above should work. You will get some 
warnings from the compiler, but the code will compile and create an executable ddscat. 

If you wish to use a different compiler (or compiler options) you will need to edit the file Make file 
to change the choice of compiler (variable FC), compilation options (variable FFLAGS), and possibly 
and loader options (variable LDFLAGS). The file Makefile as provided includes examples for com- 
pilers other than gf ortran; you may be able to simply "comment out" the section of Makefile that 
was designed for gf ortran, and "uncomment" a section designed for another compiler (e.g., Intel® 
if ort). 

6.2 Optimization 

The performance of DDSCAT 7.1 will benefit from optimization during compilation and the user should 
enable the appropriate compiler flags. 

6.3 Single vs. Double Precision 

DDSCAT 7.1 is written to allow the user to easily generate either single- or double-precision versions 
of the executable. For most purposes, the single-precision precision version of DDSCAT 7.1 should 
work fine, but if you encounter a scattering problem where the single-precision version of DDSCAT 7.1 
seems to be having trouble converging to a solution, you may wish to try using the double-precision 
version - this can be beneficial in the event that round-off error is compromising the performance of 
the conjugate-gradient solver. Of course, the double precision version will demand about twice as much 
memory as the single-precision version, and will take somewhat longer per iteration. 
The only change required is in the Makefile: for single-precision, set 

PRECISION = sp 
or for double-precision, set 

PRECISION = dp 

After changing the PRECISION variable in the Makefile (either sp — » dp, or dp — > sp), it is necessary 
to recompile the entire code. Simply type 

make clean 

make ddscat 
to create ddscat with the appropriate precision. 

6.4 OpenMP 

OpenMP is a standard for support of shared-memory parallel programming, and can provide a perfor- 
mance advantage when using DDSCAT 7.1 on platforms with multiple cpus or multiple cores. OpenMP 
is supported by many common compilers, such as gf ortran and Intel® if ort. 

If you are using a multi-cpu (or multi-core) system with OpenMP (www . openmp . org) installed, 
you can compile DDSCAT 7.1 with OpenMP directives to parallelize some of the calculations. To do 
so, simply change 

DOMP = 

OPENMP = 

to 

DOMP = -Dopenmp 
OPENMP = -openmp 
Note: OPENMP is compiler-dependent: gf ortran, for instance, requires 

OPENMP = -fopenmp 
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6.5 MKL: the Intel® Math Kernel Library 

Intel® offers the Math Kernel Library (MKL) with the ifort compiler. This library includes DFTI for 
computing FFTs. At least on some systems, DFTI offers better performance than the GPFA package. 
To use the MKL library routine DFTI: 

• You must obtain the routine mkl.dfti . f 90 and place a copy in the directory where you are 
compiling DDSCAT 7.1. mkl.df ti . f 90 is Intel® proprietary software, so we cannot distribute 
it with the DDSCAT 7.1 source code, but it should be available on any system with a license for 
the Intel® MKL library. If you cannot find it, ask your system administrator. 

• Edit the Makefile: define variables CXFFTMKL . f and CXFFTMKL . o to: 

CXFFTMFKL.f = $ (MKL_f ) 
CXFFTMKL. f = $(MKL_o) 

• Successful linking will require that the appropriate MKL libraries be available, and that the string 
LFLAGS in the Makefile be defined so as to include these libraries. 

-Imkl_em64t -lmkl_intel_thread -lmkl_core -lguide -lpthread -lmkl_intel_lp64 
appears to work on at least one installation. The Makefile contains examples. 

• type 

make clean 
make ddscat 

• The parameter file ddscat . par should have FFTMKL as the value of CMETHD. 

6.6 MPI: Message Passing Interface 

DDSCAT 7.1 includes support for parallelization under MPI. MPI (Message Passing Interface) is a 
standard for communication between processes. More than one implementation of MP I exists (e.g., 
mpich and openmpi). MPI support within DDSCAT 7.1 is compliant with the MPI-1.2 and MPI- 
2 standards^], and should be usable under any implementation of MPI that is compatible with those 
standards. 

Many scattering calculations will require multiple orientations of the target relative to the incident 
radiation. For DDSCAT 7.1, such calculations are "embarassingly parallel", because they are carried out 
essentially independently. DDSCAT 7.1 uses MPI so that scattering calculations at a single wavelength 
but for multiple orientations can be carried out in parallel, with the information for different orientations 
gathered together for averaging etc. by the master process. 

If you intend to use DDSCAT 7.1 for only a single orientation, MPI offers no advantage for DDSCAT 
7.1 so you should compile with MPI disabled. However, if you intend to carry out calculations for 
multiple orientations, and would like to do so in parallel over more than one cpu, and you have MPI 
installed on your platform, then you will want to compile DDSCAT 7.1 with MPI enabled. 

To compile with MPI disabled: in the Makefile, set 
DMPI = 

MIP . f = mpi_fake.f90 
MPI.o = mpi_fake.o 
To compile with MPI enabled: in the Makefile, set 

DMPI = -Dmpi 
MIP . f = $ (MPI.f ) 
MPI . o = $ (MPI.o) 

and edit LFLAGS as needed to make sure that you link to the appropriate MPI library (if in doubt, consult 
your systems administrator). The Makefile in the distribution includes some examples, but library 
names and locations are often system-dependent. Please do not direct questions regarding LFLAGS to 
the authors - ask your sys-admin or other experts familiar with your installation. 

Note that the MPI-capable executable can also be used for ordinary serial calculations using a single 
cpu. 



^http : / /www . mpi-f orum. org/| 
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6.7 Device Numbers IDVOUT and IDVERR 

So far as we know, there are only one operating-system-dependent aspect of DDSCAT 7.1: the device 
number to use for "standard output". 

The variables IDVOUT and IDVERR specify device numbers for "running output" and "error mes- 
sages", respectively. Normally these would both be set to the device number for "standard output" (e.g., 
writing to the screen if running interactively). Variables IDVERR are set by DATA statements in the 
"main" program DDSCAT . f 90 and in the output routine WRIMSG (file wrimsg . f 90). The executable 
statement IDVOUT=0 initializes IDVOUT to 0. In the as-distributed version of DDSCAT . f 90, the state- 
ment 

OPEN (UNIT=IDVOUT, FILE=CFLLOG) 
causes the output to UNIT=IDVOUT to be buffered and written to the file ddscat . log_nnn, where 
nnn=0 for the first cpu, 001 for the second cpu, etc. If it is desirable to have this output unbuffered 
for debugging purposes, (so that the output will contain up-to-the-moment information) simply comment 
out this OPEN statement. 

7 Moving the Executable 

Now reposition yourself into the directory DDA (e.g., type cd . .), and copy the executable from 
src/ddscat to the DDA directory by typing 

cp src/ddscat ddscat 
This should copy the file DDA/ src/ddscat to DDA/ddscat (you could, alternatively, create a sym- 
bolic link). Similarly, copy the sample parameter file ddscat . par and the file diel . tab to the DDA 
directory by typing 

cp doc/ddscat . par ddscat. par 

cp doc/diel.tab diel. tab 

8 The Parameter File ddscat. par 

When the tarfile is unpacked, it will create four directories: src, doc, diel, and examples.exp. 

The examples.exp directory has a number of subdirectories, each with files for a sample calculation. 

Here we focus on the files in the subdirectory RCTGLPRSM. 

The file ddscat . par (see also Appendix [Ab provides parameters to the program ddscat: 

' ========= Parameter file for v7 . 1 ===================' 

' -k-k-k* Preliminaries ****' 

' NOTORQ' = CMTORQ*6 (NOTORQ, DOTORQ) — either do or skip torque calculations 
'PBCGS2' = CMDSOL*6 (PBCGS2, PBCGST, PETRKP) — select solution method 
'GPFAFT' = CMDFFT* 6 (GPFAFT, FFTMKL) — select FFT method 

'GKDLDR' = CALPHA*6 (GKDLDR, LATTDR) — select prescription for polarizabilities 
'NOTBIN' = CBINFLAG (NOTBIN, ORIBIN, ALLBIN) — specify binary output 
r kkk-k Initial Memory Allocation kkkk' 

100 100 100 = dimensioning allowance for target generation 
' kkkk Target Geometry and Composition kkkk' 
'RCTGLPRSM' = CSHAPE*9 shape directive 
32 64 64 = shape parameters 1-3 

1 = NCOMP = number of dielectric materials 

' . . /diel/Au_evap' = file with refractive index 1 
' kkkk Error Tolerance kkkk' 

1.00e-5 = TOL = MAX ALLOWED (NORM OF | G>=AC I E>-ACA | X> ) / (NORM OF AC|E>) 

'kkkk Interaction cutoff parameter for PBC calculations kkkk' 

1.00e-2 = GAMMA (le-2 is normal, 3e-3 for greater accuracy) 

'kkkk Angular resolution for calculation of <cos>, etc. kkkk' 

0.5 = ETASCA (number of angles is proportional to [ ( 3+x) /ETASCA] " 2 ) 

'kkkk Wavelengths (micron) kkkk' 

0.5000 0.5000 1 'LIN' = wavelengths ( first , last , how many , how=LIN, INV, LOG) 
' kkkk Effective Radii (micron) kkkk ' 

0.49237 0.49237 1 'LIN' = aeff (first, last, how many, how=LIN, INV, LOG) 
'kkkk Define Incident Polarizations kkkk' 
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(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = IORTH (=1 to do only pol . state eOl; = 2 to also do orth. pol . state) 

'**** Specify which output files to write 

1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient. 

1 = IWRPOL (=0 to suppress, =1 to write ".pol" file for each (BETA, THETA) 
' Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 

0. 0. 1 = THETMI, THETMX, NTHETA (theta=angle between al and k) 

0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 

'**** Specify first I WAV, IRAD, IORI (normally 0) ****' 

= first I WAV, first IRAD, first IORI (0 to begin fresh) 

' **** Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 
11 12 21 22 31 41 = indices ij of elements to print 
' **** Specify Scattered Directions ****' 

'LFRAME' = CMDFRM (LFRAME, T FRAME for Lab Frame or Target Frame) 

2 = NPLANES = number of scattering planes 

0. 0. 180. 5 = phi, thetan_min, thetan_max, dtheta (in deg) for plane 1 
90. 0. 180. 5 = phi, thetan_min, thetan_max, dtheta (in deg) for plane 2 

Here we discuss the general structure of ddscat . par (see also Appendix lAb. 

8.1 Preliminaries 

ddscat . par starts by setting the values of five strings: 

• CMDTRQ specifying whether or not radiative torques are to be calculated (e.g., NOTORQ) 

• CMDSOL specifying the solution method (e.g., PBCGS2) (see section §T\\ . 

• CMDFFT specifying the FFT method (e.g., GPFAFT) (see section fl2l 

• CALPHA specifying the polarizability prescription (e.g., GKDLDR) (see TT"3l l. 

• CBINFLAG specifying whether to write out binary files (e.g., NOTBIN) 



8.2 Initial Memory Allocation 

Three integers 

MXNX MXNY MXNZ 

are given that need to be equal to or larger than the anticipated target size (in lattice units) in the xtf, 
Ytf, and ztf directions. This is required for initial memory allocation. A value like 
100 100 100 

would require initial memory allocation of only ~ 100 MBytes. Note that after the target geometry has 
been determined, DDSCAT 7.1 will proceed to reallocate only as much memory as is actually required 
for the calculation. 

8.3 Target Geometry and Composition 

• CSHAPE specifies the target geometry (e.g., RCTGLPRSM) 

As provided, the file ddscat . par is set up to calculate scattering by a 32x64x64 rectangular array of 
131072 dipoles. 

The user must specify NCOMP, which specifies how many different dielectric functions will be used. 
This is then followed by NCOMP lines, with each line giving the name (in quotes) of each dielectric 
function file. In our example, NCOMP is set to 1. The dielectric function of the target material is provided 
in the file . . /diel/Au_evap, which is gives the refractive index of evaporated Au over a range of 
wavelengths. 

8.4 Error Tolerance 

TOL = the error tolerance. Conjugate gradient iteration will proceed until the linear equations are solved 
to a fractional error TOL. The sample calculation has TOL=10 -5 . 
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8.5 Interaction Cutoff Parameter for PBC Calculations 

GAMMA = parameter limiting certain summations that are required for periodic targets (see Draine & 
Flatau 2008). The value of GAMMA has no effect on computations for finite targets - it can be set to 
any value, including 0). 

For targets that are periodic in 1 or 2 dimensions, GAMMA needs to be small for high accuracy, but 
the required cpu time increases as GAMMA becomes smaller. GAMMA = 10~ 2 is reasonable for initial 
calculations, but you may want to experiment to see if the results you are interested in are sensitive to 
the value of GAMMA. Draine & Flatau (2008) show examples of how computed results for scattering can 
depend on the value of 7. 

8.6 Angular Resolution for Computation of (cos 6), etc. 

The parameter ETASCA determines the selection of scattering angles used for computation of certain 
angular averages, such as (cos 9) and (cos 2 9), and the radiation pressure force (see Cl5l l and radia- 
tive torque (if CMDTRQ=DOTORQ). Small values of ETASCA result in increased accuracy but also cost 
additional computation time. ETASCA=0.5 generally gives accurate results. 

If accurate computation of (cos 9) or the radiation pressure force is not required, the user can set 
ETASCA to some large number, e.g. 10, to minimize unnecessary computation. 

8.7 Wavelength 

Wavelengths A are specified one line in ddscat . par consisting of values for 4 variables: 

WAVINI WAVEND NWAV CDIVID 
where WAVINI and WAVEND are real numbers, NWAV is an integer, and CDIVID is a character variable. 

• If CDIVID = ' LIN' , the A will be uniformly spaced beween WAVINI and WAVEND. 

• If CDIVID = ' INV , the A will be uniformly spaced in 1/A between WAVINI and WAVEND 

• If CDIVID = ' LOG' , the A will be uniformly spaced in log(A) between WAVINI and WAVEND 

• If CDIVID = ' TAB' , the A will be read from a user-supplied file wave . tab, with one wave- 
length per line. For this case, the values of WAVINI, WAVEND, and NWAV will be disregarded. 

The sample ddscat . par file specifies that the calculations be done for a single (A = 0.50). The units 
must be the same as the wavelength units used in the file specifying the refractive index. In this case, we 
are using fim. 

8.8 Target Size a e g 

Note that in DDSCAT the "effective radius" a e ff is the radius of a sphere of equal volume - i.e., a sphere 
of volume Nd 3 , where d is the lattice spacing and N is the number of occupied (i.e., non-vacuum) lattice 
sites in the target. Thus the effective radius a e ff = (3iV/47r) 1 / 3 (i . Our target should have a thickness 
a = 0.5/im in the xtf direction. If the rectangular solid is a x b x c, with a : b : c :: 32 : 64 : 64, then 
V = abc = 4a 3 . Thus a cS = (W/Air) 1 / 3 = (3/ 7 r) 1 /3 a = 0.49237/jm. 

8.9 Incident Polarization 

The incident radiation is always assumed to propagate along the x axis in the "Lab Frame". The sample 
ddscat .par file specifies incident polarization state e i to be along the y axis (and consequently 
polarization state 602 will automatically be taken to be along the z axis). IORTH=2 in ddscat .par 
calls for calculations to be carried out for both incident polarization states (eoi and eo2 - see Sj22]i. 
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8.10 Output Files to Write 

The sample ddscat .par file sets 1 = IWRKSC, so that ".sea" files will be written out with scattering 
information for each orientation. 

The sample ddscat .par file sets l = IWRPOL, so that a ".pol" file will be writtten out for each 
incident polarization and each orientation. Note that each such file has a size that is proportional to the 
number of dipoles, and can be quite large. For this modest test problem, with N = 131072 dipoles, the 
.pol files are each 11.8 MBytes. These files are only useful if post-processing of the solution is desired, 
such as evaluation of the electromagnetic field near the target using the DDfield program as described 
in i27] 

8.11 Target Orientation 

The target is assumed to have two vectors ai and a.2 embedded in it; a.2 is perpendicular to ai . For the 
present target shape RCTGLPRSM, the vector ai is along the xtf axis of the target, and the vector a2 is 
along the yxF axis (see ^19.1 8b . The target orientation in the Lab Frame is set by three angles: /?, 9, and 
$, defined and discussed below in CP71 Briefly, the polar angles 9 and $ specify the direction of ai in 
the Lab Frame. The target is assumed to be rotated around ai by an angle /?. The sample ddscat . par 
file specifies (3 = and $ = (see lines in ddscat .par specifying variables BETA and PHI), and 
calls for three values of the angle 9 (see line in ddscat .par specifying variable THETA). DDSCAT 
chooses 9 values uniformly spaced in cos 9; thus, asking for three values of 9 between and 90° yields 
9 = 0, 60°, and 90°. 

8.12 Starting Values of I WAV, I RAD, IORI 

Normally we begin the calculation with IWAVE=0, IRAD=0, and I OR 1=0. However, under some cir- 
cumstances, a prior calculation may have completed some of the cases. If so, the user can specify starting 
values of I WAV, IRAD, IORI; the computations will begin with this case, and continue. 

8.13 Which Mueller Matrix Elements? 

The sample parameter file specifies that DDSCAT should calculate 6 distinct Mueller scattering matrix 
elements Sy, with 11, 12, 21, 22, 31,41 being the chosen values of ij. 

8.14 What Scattering Directions? 

8.14.1 Isolated Finite Targets 

For finite targets, the user may specify the scattering directions in either the Lab Frame ( ' L FRAME ' ) or 
the Target Frame ( ' T FRAME ' ). 

For finite targets, such as specified in this sample ddscat .par, the <Sy are to be calculated for 
scattering directions specified by angles (6, <j>). 

The sample ddscat .par specifies that 2 scattering planes are to be specified: the first has = 
and the second has <fr = 90°; for each scattering plane 8 values run from to 180° in increments of 10°. 

8.14.2 1-D Periodic Targets 

For periodic targets, the scattering directions must be specified in the Target Frame: ' T FRAME ' = 
CMDFRM . 

Scattering from 1-d periodic targets is discussed in detail by Draine & Flatau (2008). For periodic 
targets, the user does not specify scattering planes. For 1 -dimensional targets, the user specifies scattering 
cones, corresponding to different scattering orders M. For each scattering cone M, the user specifies 
Cmin, Cmax> A£ ( m degrees); the azimuthal angle £ will run from £ m j n to Cmax, in increments of A£. For 
example: 
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' TFRAME ' = CMDFRM (LFRAME, T FRAME for Lab Frame or Target Frame) 
1 = number of scattering cones 

0. 0. 180. 0.05 = OrderM zetamin zetamax dzeta for scattering cone 1 

8.14.3 2-D Periodic Targets 

For targets that are periodic in 2 dimensions, the scattering directions must be specified in the Target 
Frame: ' TFRAME' = CMDFRM. 

Scattering from targets that are periodic in 2 dimensions is discussed in detail by Draine & Flatau 
(2008). For 2-D periodic targets, the user specifies the diffraction orders (M, N) for transmitted radi- 
ation: the code will automatically calculate the scattering matrix elements (M, N) for both trans- 
mitted and reflected radiation for each (M, N) specified by the user. For example: 

'TFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame) 

1 = number of scattering orders 

. . = OrderM OrderN for scattered radiation 

9 Running DDSCAT 7.1 Using the Sample dds cat . par File 

9.1 Single-Process Execution 

To execute the program on a UNIX system (running either sh or csh), simply type 

. /ddscat >& ddscat.out & 
which will redirect the "standard output" to the file ddscat . out, and run the calculation in the back- 
ground. The sample calculation [32x24x16=12288 dipole target, 3 orientations, two incident polar- 
izations, with scattering (Mueller matrix elements calculated for 38 distinct scattering directions], 
requires 27.9 cpu seconds to complete on a 2.0 GHz single-core Xeon workstation. 

9.2 Code Execution Under MPI 

Local installations of MPI will vary - you should consult with someone familiar with way MPI is in- 
stalled and used on your system. 

At Princeton University Dept. of Astrophysical Sciences we use PBS (Portable Batch Systemf|to 
schedule jobs. MPI jobs are submitted using PBS by first creating a shell script such as the following 
example file pbs . submit: 

#!/bin/bash 

#PBS -1 nodes=2:ppn=l 

#PBS -1 mem=1200MB,pmem=300MB 

#PBS -m bea 

#PBS -j oe 

cd $PBS_0_WORKDIR 
/usr/local/bin/mpiexec ddscat 

The lines beginning with #PBS -1 specify the required resources: 

#PBS -1 nodes=2 :ppn=l specifies that 2 nodes are to be used, with 1 processor per node. 
#PBS -1 mem=1200MB, pmem=300MB specifies that the total memory required (mem) is 1200MB, 
and the maximum physical memory used by any single process (pmem) is 300MB. The actual definition 
of mem is not clear, but in practice it seems that it should be set equal to 2x(nodes)x(ppn)x(pmem) . 
#PBS -m bea specifies that PBS should send email when the job begins (b), and when it ends (e) or 
aborts (a). 

#PBS -j oe specifies that the output from stdout and stderr will be merged, intermixed, as stdout. 



" jhttp : / / www . openpbs . org| 
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This example assumes that the executable dscat is located in the same directory where the code is to 
execute and write its output. If ddscat is located in another directory, simply give the full pathname, 
to it. The qsub command is used to submit the PBS job: 

qsub pbs. submit 

As the calculation proceeds, the usual output files will be written to this directory: for each wave- 
length, target size, and target orientation, there will be a file waaarbbb~kccc .sea, where aaa=000, 001, 
002, ... specifies the wavelength, bbb=000, 001, 002, ... specifies the target size, and cc=000, 001, 002, ... 
specifies the orientation. For each wavelength and target size there will also be a file waaarbbbori . avg 
with orientationally-averaged quantities. Finally, there will also be tables qtable, and qtable2 with 
orientationally-averaged cross sections for each wavelength and target size. 

In addition, each processor employed will write to its own log file ddscat . log jinn, where nnn=000, 
001, 002, ... . These files contain information concerning cpu time consumed by different parts of the cal- 
culation, convergence to the specified error tolerance, etc. If you are uncertain about how the calculation 
proceeded, examination of these log files is recommended. 

10 Output Files 
10.1 ASCII files 

If you run DDSCAT using the command 
ddscat >& ddscat. out & 
you will have various types of ASCII files when the computation is complete: 

• a file ddscat . out; 

• a file ddscat . log_000; 

• a file mt able; 

• a file qtable; 

• a file qtable2; 

• files wxxxryyyori . avg (one, wOOOrOOOori . avg, for the sample calculation); 

• if ddscat . par specified IWRKSC=1, there will also be files wxxxryyykzzz . sea (3 for the sam- 
ple calculation: wOOOrOOOkOOO . sea, wOOOrOOOkOOl . sea, w000r000k002 . sea). 

• if ddscat . par specified IWRPOL=l, there will also be files wxxxryyykzzz . pol« for n=l and 
(if IORTH=2)«=2 

The file ddscat . out will contain minimal information (it may in fact be empty). 

The file ddscat . log_000 will contain any error messages generated as well as a running report 
on the progress of the calculation, including creation of the target dipole array. During the iterative 
calculations, Q C xt, Qabs, and Q p h a are printed after each iteration; you will be able to judge the degree 
to which convergence has been achieved. Unless TIME IT has been disabled, there will also be timing 
information. If the MP I option is used to run the code on multiple cpus, there will be one file of the form 
ddscat . log_nnn for each of the cpus, with nnn=000 , 001,002, . . .. 

The file mtable contains a summary of the dielectric constant used in the calculations. 

The file qtable contains a summary of the orientationally-averaged values of Q C xt, Qabs, Qsca, 
g(l) = (cos(6> s )), (cos 2 (0 s )), Qbk, and iV sca . Here Q oxt , Qabs, and Q sca are the extinction, absorption, 
and scattering cross sections divided by 7ra;f ff . Qbk is the differential cross section for backscattering 
(area per sr) divided by tto^ s . N sca is the number of scattering directions used for averaging over 
scattering directions (to obtain (cos 0), etc.) (see < j23l . 

The file qtable2 contains a summary of the orientationally-averaged values of Q p h a , Qpoh and 
Qcpoi- Here Q p i ia is the "phase shift" cross section divided by 7ra^ ff (see definition in Draine 1988). 
Qpoi is the "polarization efficiency factor", equal to the difference between Q cx t for the two orthogonal 
polarization states. We define a "circular polarization efficiency factor" Q cpo i = Q po iQpha, since an 
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optically-thin medium with a small twist in the alignment direction will produce circular polarization in 
initially unpolarized light in proportion to Q cpo i- 

For each wavelength and size, DDSCAT 7.1 produces a file with a name of the form 
vrxxxryyyori . avg, where index xxx (=000, 001, 002....) designates the wavelength and index yyy 
(=000, 001, 002...) designates the "radius"; this file contains Q values and scattering information aver- 
aged over however many target orientations have been specified (see CP71 The file wOOOrOOOori.avg 
produced by the sample calculation is provided below in AppendixlBl 

In addition, if ddscat .par has specified IWRKSC=1 (as for the sample calculation), DDSCAT 
7.1 will generate files with names of the form wxxxryyykzzz ■ avg, where xxx and yyy are as before, and 
index zzz =(000,001,002...) designates the target orientation; these files contain Q values and scattering 
information for each of the target orientations. The structure of each of these files is very similar to 
that of the wxxxryyyor ± . avg files. Because these files may not be of particular interest, and take up 
disk space, you may choose to set IWRKSC=0 in future work. However, it is suggested that you run the 
sample calculation with IWRKSC=1. 

The sample ddscat . par file specifies IWRKSC=1 and calls for use of 1 wavelength, 1 target size, 
and averaging over 3 target orientations. Running DDSCAT 7.1 with the sample ddscat . par file will 
therefore generate files wOOOrOOOkOOO. sea, wOOOrOOOkOOl. sea, and w000r000k002 . sea . 
To understand the information contained in one of these files, please consult AppendixICl which contains 
an example of the file wOOOrOOOkOOO.sca produced in the sample calculation. 

10.2 Binary Option 

It is possible to output an "unformatted" or "binary" file (dd.bin) with fairly complete information, 
including header and data sections. This is accomplished by specifying either ALLBIN or ORIBIN in 
ddscat . par . 

Subroutine wr itebin . f 9 provides an example of how this can be done. The "header" section 
contains dimensioning and other variables which do not change with wavelength, particle geometry, 
and target orientation. The header section contains data defining the particle shape, wavelengths, parti- 
cle sizes, and target orientations. If ALLBIN has been specified, the "data" section contains, for each 
orientation, Mueller matrix results for each scattering direction. The data output is limited to actual di- 
mensions of arrays; e.g. nscat ,4,4 elements of Muller matrix are written rather than mxscat ,4,4. 
This is an important consideration when writing postprocessing codes. 



11 Choice of Iterative Algorithm 

As discussed elsewhere (e.g., Draine 1988), the problem of electromagnetic scattering of an incident 
wave Einc by an array of N point dipoles can be cast in the form 

AP = E (17) 

where E is a 3iV-dimensional (complex) vector of the incident electric field Ej nc at the N lattice sites, 
P is a 3./V-dimensional (complex) vector of the (unknown) dipole polarizations, and A is a 3iV x 37V 
complex matrix. 

Because 37V is a large number, direct methods for solving this system of equations for the unknown 
vector P are impractical, but iterative methods are useful: we begin with a guess (typically, P = 0) for 
the unknown polarization vector, and then iteratively improve the estimate for P until equation (fTTI i is 
solved to some error criterion. The error tolerance may be specified as 

|AtAP-A+E| , 

[ < h , (18) 



|AtE| 

where A^ is the Hermitian conjugate of A [(A')ij = (A/;)*], and h is the error tolerance. We typically 
use h = 10~ 5 in order to satisfy eq.JTTb to high accuracy. The error tolerance h can be specified by the 
user through the parameter TOL in the parameter file ddscat . par (see AppendixlAt. 
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A major change in going from DDSCAT.4b to 5a (and subsequent versions) was the implementation 
of several different algorithms for iterative solution of the system of complex linear equations. DDSCAT 
7.1 is now structured to permit solution algorithms to be treated in a fairly "modular" fashion, facilitating 
the testing of different algorithms. A number of algorithms were compared by Flatau (1997^3; two of 
them (PBCGST and PETRKP) performed well and are made available to the user in DDSCAT. DDSCAT 
7.1 now also offers a third option, PBCGS2. The choice of algorithm is made by specifying one of the 
options: 

• PBCGS2 - BiConjugate Gradient with Stabilization as implemented in the routine ZBCG2 by 
M.A. Botchev, University of Twente. This is based on the PhD thesis of D.R. Fokkema, and on 
work by Sleijpen & vanderVorst (1995,1996). 

• PBCGST - Preconditioned BiConjugate Gradient with STabilitization method from the Parallel 
Iterative Methods (PIM) package created by R. Dias da Cunha and T. Hopkins. 

• PETRKP - the complex conjugate gradient algorithm of Petravic & Kuo-Petravic (1979), as coded 
in the Complex Conjugate Gradient package (CCGPACK) created by P.J. Flatau. This is the algo- 
rithm discussed by Draine (1988) and used in previous versions of DDSCAT. 

All three methods work fairly well. Our experience suggests that PBCGS2 is generally fastest and best- 
behaved, and we recommend that the user try it first. However, in case it runs into numerical difficulties, 
PBCGST and PETRKP are available as alternatives]! 



12 Choice of FFT Algorithm 

DDSCAT 7.1 offers two FFT options: (1) the GPFA FFT algorithm developed by Dr. Clive Temperton 
(1992),0 and (2) the Intel® MKL routine DFTI. 

A postscript copy of this report - file eg . ps - is distributed with the DDSCAT 7.1 documentation. 

6 The Parallel Iterative Methods (PIM) by Rudnei Dias da Cunha (rddSukc . ac . uk) and Tim Hopkins (trh@ukc . ac . uk) 
is a collection of Fortran 77 routines designed to solve systems of linear equations on parallel and scalar computers using a variety 
of iterative methods (available at 

|http : / /www . mat ■ uf rgs . br /pim-e . html) . PIM offers a number of iterative methods, including 

• Conjugate-Gradients, CG (Hestenes 1952), 

• Bi-Conjugate-Gradients, BICG (Fletcher 1976), 

• Conjugate-Gradients squared, CGS (Sonneveld 1989), 

• the stabilised version of Bi-Conjugate-Gradients, BICGSTAB (van der Vorst 1991), 

• the restarted version of BICGSTAB, RBICGSTAB (Sleijpen & Fokkema 1992) 

• the restarted generalized minimal residual, RGMRES (Saad 1986), 

• the restarted generalized conjugate residual, RGCR (Eisenstat 1983), 

• the normal equation solvers, CGNR (Hestenes 1952 and CGNE (Craig 1955), 

• the quasi-minimal residual, QMR (highly parallel version due to Bucker & Sauren 1996), 

• transpose-free quasi-minimal residual, TFQMR (Freund 1992), 

• the Chebyshev acceleration, CHEBYSHEV (Young 1981). 

The source code for these methods is distributed with DDSCAT but only PBCGST and PETRKP can be called directly via 
ddscat . par. It is possible to add other options by changing the code in get f ml. f 90 . Flatau (1997) has compared 
the convergence rates of a number of different methods. A helpful introduction to conjugate gradient methods is provided 
by the report "Conjugate Gradient Method Without Agonizing Pain" by Jonathan R. Shewchuk, available as a postscript file: 
|ftp : //REPORTS .ADM. CS . CMU . EDU/usr / anon/ 1 9 94 /CMU-CS-94-12 5 . ps) 

7 The GPFA code contains a parameter LVR which is set in data statements in the routines gpf a2f , gpf a3f , and gpf a5f . 
LVR is supposed to be optimized to correspond to the "length of a vector register" on vector machines. As delivered, this parameter 
is set to 64, which is supposed to be appropriate for Crays other than the C90. For the C90, 128 is supposed to be preferable (and 
perhaps "preferable" should be read as "necessary" - there is some basis for fearing that results computed on a C90 with LVR other 
than 128 run the risk of being incorrect!) The value of LVR is not critical for scalar machines, as long as it is fairly large. We found 
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Figure 5: Comparison of cpu time required by 3 different FFT implementations. It is seen that the GPFA and FFTW implemen- 
tations have comparable speeds, much faster than Brenner's FFT implementation. 

The GPFA routine is portable and quite fast: Figure 5 compares the speed of three FFT implemen- 
tations: Brenner's, GPFA, and FFTW Qhttp : / /www . f f tw . org) . We see that while for some cases 
FFTW 2 . 1 . 5 is faster than the GPFA algorithm, the difference is only marginal. The FFTW code and 
GPFA code are quite comparable in performance - for some cases the GPFA code is faster, for other 
cases the FFTW code is faster. For target dimensions which are factorizable as 2 8 3 ;? 5 (for integer i, 
j, k), the GPFA and FFTW codes have the same memory requirements. For targets with extents N x , 
N y , N z which are not factorizable as 2 i 3 : '5 fe , the GPFA code needs to "extend" the computational vol- 
ume to have values of N x , N y , and N z which are factorizable by 2, 3, and 5. For these cases, GPFA 
requires somewhat more memory than FFTW. However, the fractional difference in required memory is 
not large, since integers factorizable as 2 I 3 : '5 fc occur fairly frequently!! [Note: This "extension" of the 
target volume occurs automatically and is transparent to the user.] 

DDSCAT 7.1 offers a new FFT option: the Intel® Math Kernel Library DFTI. This is tuned for 
optimum performance, and appears to offer real performance advantages on modern multi-core cpus. 
With this now available, the FFTW option, which had been included in DDSCAT 6.1, has been removed 
from DDSCAT 7.1. 

The choice of FFT implementation is obtained by specifying one of: 

• FFTMKL to use the Intel® MKL routine DFTI (see %.5i . This is recommended, but requires 
that the Intel® Math Kernel Library be installed on your system. 

little difference between LVR=64 and 128 on a Sparc 10/51, on an Ultrasparc 170, and on an Intel® Xeon cpu. You may wish to 
experiment with different LVR values on your computer architecture. To change LVR, you need to edit gpf a . f 90 and change the 
three data statements where LVR is set. 

8 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 
125, 128, 135, 144, 150, 160, 162, 180, 192, 200 216, 225, 240, 243, 250, 256, 270, 288, 300, 320, 324, 360, 375, 384, 400, 405, 
432, 450, 480, 486, 500, 512, 540, 576, 600, 625, 640, 648, 675, 720, 729, 750, 768, 800, 810, 864, 900, 960, 972, 1000, 1024, 
1080, 1125, 1152, 1200, 1215, 1250, 1280, 1296, 1350, 1440, 1458, 1500, 1536, 1600, 1620, 1728, 1800, 1875, 1920, 1944, 2000, 
2025, 2048, 2160, 2187, 2250, 2304, 2400, 2430, 2500, 2560, 2592, 2700, 2880, 2916, 3000, 3072, 3125, 3200, 3240, 3375, 3456, 
3600, 3645, 3750, 3840, 3888, 4000, 4050, 4096 are the integers < 4096 which are of the form 2 i 3 j 5 fc . 
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Figure 6: Scattering and absorption for am — 1.7 + O.li sphere, calculated using two prescriptions for the polarizability: 
LATTDR is the lattice dispersion relation result of Draine & Goodman (1993). GKDLDR is the lattice dispersion relation result 
of Gutkowicz-Krusin & Draine (2004). Results are shown as a function of scattering parameter x — 2ira/\; the upper scale gives 
values of |m|fcd. We see that the cross sections calculated with these two prescriptions are quite similar for |m|fcd ^, 0.5. For other 
examples see Gutkowicz-Krusin & Draine (2004). 

• GPFAFT to use the GPFA algorithm (Temperton 1992).This is is not quite as fast as FFTMKL, 
but is written in plain Fortran-90. 

13 Dipole Polarizabilities 

Option GKDLDR specifies that the polarizability be prescribed by the "Lattice Dispersion Relation", 
with the polarizability found by Gutkowicz-Krusin & Draine (2004), who corrected a subtle error in the 
analysis of Draine & Goodman (1993). For \m\kd £ 1, the GKDLDR polarizability differs slightly from 
the LATTDR polarizability, but the differences in calculated scattering cross sections are relatively small, 
as can be seen from Figure[6] We recommend option GKDLDR. 

Users wishing to compare can invoke option LATTDR to specify that the "Lattice Dispersion Rela- 
tion" of Draine & Goodman (1993) be employed to determine the dipole polarizabilities. This polariz- 
ability also works well. 

14 Dielectric Functions 

In order to assign the appropriate dipole polarizabilities, DDSCAT 7.1 must be given the dielectric 
constant of the material (or materials) of which the target of interest is composed relative to the dielectric 
constant of the ambient medium. This information is supplied to DDSCAT 7.1 through a table (or tables), 
read by subroutine DIELEC in file dielec . f 90, and providing either the complex refractive index 
m = n + ik or complex dielectric function e = ei + ie 2 as a function of wavelength A. Since m = e 1 / 2 , 
or e = m 2 , the user must supply either more. 

DDSCAT 7.1 can calculate scattering and absorption by targets with anisotropic dielectric functions, 
with arbitrary orientation of the optical axes relative to the target shape. See £1261 

As discussed in £13.21 DDSCAT calculates scattering for targets embedded in dielectric media us- 
ing the relative dielectric constant e = etargct(w)/e me dium(w) or the relative refractive index m = 
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WtargetOj-O/mmediumO^) and the wavelength in the medium A = (27rc/w)/m me dium> where c is the 
speed of light in vacuo. 

It is therefore important that the table containing the dielectric function information give the relative 
dielectric function or refractive index as a function of the wavelength in the ambient medium. 

The table formatting is intended to be quite flexible. The first line of the table consists of text, up 
to 80 characters of which will be read and included in the output to identify the choice of dielectric 
function. (For the sample problem, it consists of simply the statement m = 1.33 + O.Oli.) The 
second line consists of 5 integers; either the second and third or the fourth and fifth should be zero. 

• The first integer specifies which column the wavelength is stored in. 

• The second integer specifies which column Re(m) is stored in. 

• The third integer specifies which column Im(m) is stored in. 

• The fourth integer specifies which column Re(e) is stored in. 

• The fifth integer specifies which column Im(e) is stored in. 

If the second and third integers are zeros, then DIELEC will read Re(e) and Im(e) from the file; if the 
fourth and fifth integers are zeros, then Re(m) and Im(m) will be read from the file. 

The third line of the file is used for column headers, and the data begins in line 4. There must be at 
least 3 lines of data: even if e is required at only one wavelength, please supply two additional "dummy" 
wavelength entries in the table so that the interpolation apparatus will not be confused. 



15 Calculation of (cos 9), Radiative Force, and Radiation Torque 

In addition to solving the scattering problem for a dipole array, DDSCAT can compute the three- 
dimensional force F ra d and torque r rac j exerted on this array by the incident and scattered radiation 
fields. The radiation torque calculation is carried out, after solving the scattering problem, only if 
DOTORQ has been specified in ddscat .par. For each incident polarization mode, the results are 
given in terms of dimensionless efficiency vectors Q pr and Qr, defined by 

Qpr = 2 rad , (19) 
™ cff Urad 

Qr ^ J£^L_ (20 ) 

7ra; ff w ra d 

where F ra d and r ra( j are the time-averaged force and torque on the dipole array, k = 2ir/\ is the 
wavenumber in vacuo, and u Ta d = Eq/8it is the time-averaged energy density for an incident plane 
wave with amplitude Eq cos(wi + cf>). The radiation pressure efficiency vector can be written 

Qpr = Qoxt.k - Qscag , (21) 

where k is the direction of propagation of the incident radiation, and the vector g is the mean direction 
of propagation of the scattered radiation: 

1 f , n dC sca (n,k) „ 

g= c^r ^— n ' (22) 

where dfl is the element of solid angle in scattering direction n, and dC sca /d£l is the differential scat- 
tering cross section. The components of Q pr are reported in the Target Frame: Q pr ,i = F ra d • xtf, 

Qpr, 2 = F ra d • YTF, Qpr, 3 = F ra d ■ ZTF- 

Equations for the evaluation of the radiative force and torque are derived by Draine & Weingartner 
(1996). It is important to note that evaluation of Q pr and Qr involves averaging over scattering directions 
to evaluate the linear and angular momentum transport by the scattered wave. This evaluation requires 
appropriate choices of the parameter ETASCA - see S|23] 

In addition, DDSCAT calculates (cos 9) [the first component of the vector g in eq. d22b l and the 
second moment (cos 2 9). These two moments are useful measures of the anisotropy of the scattering. 
For example, Draine (2003) gives an analytic approximation to the scattering phase function of dust 
mixtures that is parameterized by the two moments (cos 9) and (cos 2 9). 
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16 Memory Requirements 

The memory requirements are determined by the size of the "computational volume" - this is a rectan- 
gular region, of size NXxNYxNZ that is large enough to contain the target. If using the GPFAFT option, 
then NX, NY, NZ are also required to have only 2,3, and 5 as prime factors (see footnote[8]). 
In single precision, the memory requirement for DDSCAT 7.1 is approximately 

(35. + 0.0010 x NX x NY x NZ) Mbytes for single precision (23) 

(42 + 0.0020 x NX x NY x NZ) Mbytes for double precision (24) 

Thus, in single precision, a 48x48x48 calculation requires ~146 MBytes. 

The memory is allocated dynamically - once the target has been created, DDSCAT 7.1 will deter- 
mine just how much overall memory is needed, and will allocate it. However, the user must provide 
information (via dels cat . par) to allow DDSCAT 7.1 to allocate sufficiently large arrays to carry out 
the initial target creation. Initially, the only arrays that will be allocated are those related to the target 
geometry, so it is OK to be quite generous in this initial allowance, as the memory required for the target 
generation step is small compared to the memory required to carry out the full scattering calculation. 

17 Target Orientation 

Recall that we define a "Lab Frame" (LF) in which the incident radiation propagates in the +x direction. 
For purposes of discussion we will always let unit vectors xlf, Ylf, zlf = xlf x Ylf be the three 
coordinate axes of the LF. 

In ddscat .par one specifies the first polarization state eoi (which obviously must lie in the y, z 
plane in the LF); DDSCAT automatically constructs a second polarization state eo2 = xlf x &%\ orthog- 
onal to eoi. Users will often find it convenient to let polarization vectors eoi = y> ^02 = z (although 
this is not mandatory - see ! j22l . 

Recall that definition of a target involves specifying two unit vectors, ai and a.2, which are imagined 
to be "frozen" into the target. We require a.2 to be orthogonal to ai . Therefore we may define a "Target 
Frame" (TF) defined by the three unit vectors ai, a.2, and EL3 = ai x a.2 . 

For example, when DDSCAT creates a 32x24x16 rectangular solid, it fixes ai to be along the 
longest dimension of the solid, and a.2 to be along the next-longest dimension. 

Important Note: for periodic targets, DDSCAT 7.1 requires that the periodic target have ai = xtf 
and a 2 = ytf- 

Orientation of the target relative to the incident radiation can in principle be determined two ways: 

1 . specifying the direction of ai and a.2 in the LF, or 

2. specifying the directions of xlf (incidence direction) and yLF in the TF. 

DDSCAT uses method 1.: the angles 9, $, and /3 are specified in the file ddscat . par. The target is 
oriented such that the polar angles 9 and <E> specify the direction of ai relative to the incident direction 
xlf, where the xlfjYlf plane has $ = 0. Once the direction of ai is specified, the angle (3 then 
specifies how the target is to rotated around the axis ai to fully specify its orientation. A more extended 
and precise explanation follows: 

17.1 Orientation of the Target in the Lab Frame 

DDSCAT uses three angles, 9, <&, and /3, to specify the directions of unit vectors ai and a.2 in the LF 
(see Fig. [7}. 

9 is the angle between ai and xlf- 

When $ = 0, ai will lie in the xlf, Ylf plane. When $ is nonzero, it will refer to the rotation of ai 
around xlf: e-g-, $ = 90° puts a.! in the Xlf, zlf plane. 

When j3 = 0, a 2 will lie in the xlf , &i plane, in such a way that when 9 = and <I> = 0, a 2 is in 
the ylf direction: e.g, 9 = 90°, $ = 0, f3 = has = yLF an d &2 = — ^lf- Nonzero j3 introduces 
an additional rotation of a2 around ai: e.g., 9 = 90°, <& = 0, (3 = 90° has ai = yLF and a.2 = zlf- 
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Fi gure 7; Target orientation in the Lab Frame, x = xlf is the direction of propagation of the incident radiation, and y = yu 
is the direction of the real component (at xlf = 0, t = 0) of the first incident polarization mode. In this coordinate system, the 
orientation of target axis ai is specified by angles O and <E>. With target axis ai fixed, the orientation of target axis a.2 is then 
determined by angle fi specifying rotation of the target around an. When fi = 0, as lies in the ai,XLF plane. 
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Mathematically: 



ai = xlf cos 9 + yLF sin 9 cos $ + zlf sin 9 sin $ (25) 
a 2 = —xlf sin 9 cos j3 + yLF [cos 9 cos (3 cos <J> — sin (3 sin <£>] 

+zlf [cos 9 cos f3 sin $ + sin f3 cos $] (26) 
a 3 = xlf sin 9 sin j3 — yLF [cos 9 sin (3 cos $ + cos j3 sin $] 

—zlf [cos 9 sin /3 sin <E> — cos /3 cos $] (27) 



or, equivalently: 



xlf = a-i cos 9 — a 2 sin 9 cos f3 + a 3 sin 9 sin j3 (28) 
Ylf = ai sin 9 cos $ + a 2 [cos 9 cos f3 cos $ — sin /3 sin $] 

— a 3 [cos 9 sin (3 cos $ + cos (3 sin $] (29) 
zlf = &i sin 9 sin $ + a 2 [cos 9 cos (3 sin $ + sin (3 cos $] 

— a 3 [cos 9 sin (3 sin $ — cos /3 cos $] (30) 

17.2 Orientation of the Incident Beam in the Target Frame 

Under some circumstances, one may wish to specify the target orientation such that xlf (the direction 
of propagation of the radiation) and yLF (usually the first polarization direction) and zlf (= xlf x 
Ylf) refer to certain directions in the TF. Given the definitions of the LF and TF above, this is simply 
an exercise in coordinate transformation. For example, one might wish to have the incident radiation 
propagating along the (1,1,1) direction in the TF (example 14 below). Here we provide some selected 



examples: 








1. 


xlf 


= ai, yLF = 


a 2 , z L F = 


a 3 :6 = 0, $ + /3 = 


2. 


xlf 


= ai, yLF = 


a 3 , zlf = 


-a 2 : 9 = 0, $ + p = 90° 


3. 


xlf 


= a 2 , yLF = 


ax, zlf = 


-a 3 : 9 = 90°, P = 180°, $ = 


4. 


xlf 


= a 2 , yLF = 


a 3 , z L F = 


ax : 6 = 90°, P = 180°, $ = 90° 


5. 


xlf 


= a 3 , yLF = 


ax, zlf = 


a 2 : 9 = 90°, P = -90°, $ = 


6. 


xlf 


= a 3 , yLF = 


a 2 , z L F = 


-ax : 9 = 90°, P = -90°, $ = -90° 


7. 


xlf 


= -ax, yLF 


= a 2 , zlf 


= -a 3 : 9 = 180°, P + $ = 180° 


8. 


xlf 


= -ai, yLF 


= a 3 , z L F 


= a 2 : 9 = 180°, p + $ = 90° 


9. 


xlf 


= -a 2 , yLF 


= ax, z L F 


= a 3 : 9 = 90°, P = 0, $ = 


10. 


xlf 


= -a 2 , yLF 


= a 3 , z L F 


= -a x : 6 = 90°, P = 0, $ = -90° 


11. 


xlf 


= -a 3 , yLF 


= a 1; z L F 


= -a 2 : 9 = 90°, P = -90°, $ = 


12. 


xlf 


= -a 3 , yLF 


= a 2 , z L F 


= ax : 9 = 90°, P = -90°, $ = 90° 


13. 


xlf 


= (ai + a 2 )/\/2, yLF = 


= a 3 , z LF = (ai - a 2 )/72 : 9 = 45°, p 


14. 


xlf 


= (ai + a 2 - 


h a 3 )/ v/3, yLF = (ax - a 2 )/ ^/2, z L f = (ax + a 2 - 




9 = 


54.7356°, /3 


= 135°, $ 


= 30°. 



180°, $ = 90° 



17.3 Sampling in 6, $, and f3 

The present version, DDSCAT 7.1, chooses the angles /?, 9, and $ to sample the intervals (BETAMI , BETAMX), 
(THETMI , THETMX) , (PHIMIN, PHIMAX), where BETAMI, BETAMX, THETMI, THETMX, PHIMIN, 
PHIMAX are specified in ddscat . par . The prescription for choosing the angles is to: 

• uniformly sample in P; 

• uniformly sample in $; 
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• uniformly sample in cos 9. 

This prescription is appropriate for random orientation of the target, within the specified limits of f3, i>, 
and9. 

Note that when DDSCAT 7.1 chooses angles it handles (5 and $ differently from O. The range for j3 
is divided into NBETA intervals, and the midpoint of each interval is taken. Thus, if you take BETAMI=0, 
BETAMX=90,NBETA=2 you will get /3 = 22.5° and 67.5°. Similarly, if you take PHIMIN=0, PHIMAX=180, 
NPHI=2 you will get $ = 45° and 135°. 

Sampling in 9 is done quite differently from sampling in (3 and <£>. First, as already mentioned above, 
DDSCAT 7.1 samples uniformly in cos 9, not 9. Secondly, the sampling depends on whether NTHETA 
is even or odd. 

• If NTHETA is odd, then the values of 6 selected include the extreme values THETMI and THETMX; 
thus, THETMI=0, THETMX=90, NTHETA=3 will give you 6 = 0, 60°, 90°. 

• If NTHETA is even, then the range of cos 9 will be divided into NTHETA intervals, and the mid- 
point of each interval will be taken; thus, THETMI=0, THETMX=90, NTHETA=2 will give you 
6 = 41.41° and 75.52° [cos 9 = 0.25 and 0.75]. 

The reason for this is that if odd NTHETA is specified, then the "integration" over cos 9 is performed 
using Simpson's rule for greater accuracy. If even NTHETA is specified, then the integration over cos 9 
is performed by simply taking the average of the results for the different 9 values. 

If averaging over orientations is desired, it is recommended that the user specify an odd value of 
NTHETA so that Simpson's rule will be employed. 

18 Orientational Averaging 

DDSCAT has been constructed to facilitate the computation of orientational averages. How to go about 
this depends on the distribution of orientations which is applicable. 

18.1 Randomly-Oriented Targets 

For randomly-oriented targets, we wish to compute the orientational average of a quantity Q(f3, 9, $): 

(Q) = ^2j d P J dcosQJ d$Q(fi,Q,$) . (31) 

To compute such averages, all you need to do is edit the file ddscat .par so that DDSCAT knows 
what ranges of the angles j3, 9, and $ are of interest. For a randomly-oriented target with no symmetry, 
you would need to let (3 run from to 360°, 9 from to 180°, and <f> from to 360°. 

For targets with symmetry, on the other hand, the ranges of j3, 9, and <f> may be reduced. First 
of all, remember that averaging over <f> is relatively "inexpensive", so when in doubt average over to 
360°; most of the computational "cost" is associated with the number of different values of (/3,9) which 
are used. Consider a cube, for example, with axis a.i normal to one of the cube faces; for this cube /3 
need run only from to 90°, since the cube has fourfold symmetry for rotations around the axis ai. 
Furthermore, the angle 9 need run only from to 90°, since the orientation (f3,0,$) is indistinguishable 
from (/3, 180° - 9, 360° - $). 

For targets with symmetry, the user is encouraged to test the significance of /3,9,<& on targets with 
small numbers of dipoles (say, of the order of 100 or so) but having the desired symmetry. 

18.2 Nonrandomly-Oriented Targets 

Some special cases (where the target orientation distribution is uniform for rotations around the x axis = 
direction of propagation of the incident radiation), one may be able to use DDSCAT 7.1 with appropriate 
choices of input parameters. More generally, however, you will need to modify subroutine ORIENT to 
generate a list of NBETA values of /3, NTHETA values of 9, and NPHI values of <E>, plus two weight- 
ing arrays WGTA ( 1-NTHETA, 1-NPHI) and WGTB ( 1 -NBETA) . Here WGTA gives the weights which 
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should be attached to each (0,$) orientation, and WGTB gives the weight to be attached to each /3 ori- 
entation. Thus each orientation of the target is to be weighted by the factor WGTAxWGTB. For the case 
of random orientations, DDSCAT chooses 9 values which are uniformly spaced in cos 0, and j3 and $ 
values which are uniformly spaced, and therefore uses uniform weights 
WGTB=1./NBETA 

When NTHETA is even, DDSCAT sets 
WGTA=1./(NTHETAXNPHI) 

but when NTHETA is odd, DDSCAT uses Simpson's rule when integrating over and 
WGTA= (1/3 or 4/3 or2/3)/(NTHETAXNPHl) 
Note that the program structure of DDSCAT may not be ideally suited for certain highly oriented 
cases. If, for example, the orientation is such that for a given $ value only one value is possible 
(this situation might describe ice needles oriented with the long axis perpendicular to the vertical in the 
Earth's atmosphere, illuminated by the Sun at other than the zenith) then it is foolish to consider all the 
combinations of and $ which the present version of DDSCAT is set up to do. We hope to improve 
this in a future version of DDSCAT. 

19 Target Generation: Isolated Finite Targets 

DDSCAT contains routines to generate dipole arrays representing finite targets of various geometries, 
including spheres, ellipsoids, rectangular solids, cylinders, hexagonal prisms, tetrahedra, two touching 
ellipsoids, and three touching ellipsoids. The target type is specified by variable CSHAPE on line 9 of 
ddscat .par, up to 12 target shape parameters (SHPARi, SHPAR 2 , SHPAR 3 , ...) on line 10. The 
target geometry is most conveniently described in a coordinate system attached to the target which we 
refer to as the "Target Frame" (TF), with orthonormal unit vectors xtf, Ytf, ztf = xtf x Ytf- Once 
the target is generated, the orientation of the target in the Lab Frame is accomplished as described in CFTl 
Every target generation routine will specify 

• The "occupied" lattice sites; 

• The composition associated with each occupied lattice site; 

• Two "target axes" a.i and &2 that are used as references when specifying the target orientation; and 

• The location of the Target Frame origin of coordinates. 
Target geometries currently supported include: 

19.1 FROM_FILE = Target composed of possibly anisotropic material, defined 
by list of dipole locations and "compositions" obtained from a file 

If anisotropic, the "microcrystals" in the target are assumed to be aligned with the principal axes of 
the dielectric tensor parallel to Xtf, 9tf^ an d %tf- This option causes DDSCAT to read the target 
geometry and composition information from a file shape . dat instead of automatically generating 
one of the geometries for which DDSCAT has built-in target generation capability. The shape . dat 
file is read by routine REASHP (file reashp . f 90). The file shape . dat gives the number N of 
dipoles in the target, the components of the "target axes" ai and a 2 in the Target Frame (TF), the vector 
Xo(l — 3) determining the correspondence between the integers IXYZ and actual coordinates in the TF, 
and specifications for the location and "composition" of each dipole. The user can customize REASHP as 
needed to conform to the manner in which the target description is stored in file shape . dat . However, 
as supplied, REASHP expects the file shape . dat to have the following structure: 

• one line containing a description; the first 67 characters will be read and printed in various output 
statements 

• N = number of dipoles in target 

• a>ix a>\y a iz = x,y,z components (in TF) of ai 
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• o,2x a,2 V a,2z = x,y,z components (in TF) of a.2 



• d x /d dy/d d z /d = 1. 1. 1. = relative spacing of dipoles in xtf, Ytf, ztf directions 

• xq x xq v xqz = TF coordinates xtf / d yrp/d zrp/d corresponding to lattice site IXYZ= 

• (line containing comments) 

• dummy IXYZ (1, 1) IXYZ(1,2) IXYZ(1,3) IC0MP(1,1) IC0MP(1,2) IC0MP(1,3) 

• dummy IXYZ (2, 1) IXYZ (2, 2) IXYZ (2, 3) IC0MP(2,1) ICOMP(2,2) ICOMP(2,3) 

• dummy IXYZ (3, 1) IXYZ (3, 2) IXYZ (3, 3) IC0MP(3,1) ICOMP(3,2) ICOMP(3,3) 

• ... 

• dummy IXYZ (J, 1) IXYZ (J, 2) IXYZ (J, 3) ICOMP ( J, 1) IC0MP(J,2) IC0MP(J,3) 

• ... 

• dummy IXYZ (N, 1) IXYZ (N, 2) IXYZ (N, 3) ICOMP (N,l) ICOMP (N, 2) ICOMP (N, 3) 

where dummy is a number (integer or floating point)that might, for example, identify the dipole. This 
number will not be used in any calculations. 

If the target material at location J is isotropic, ICOMP ( J, 1 ) , ICOMP ( J, 2 ) , and ICOMP ( J, 3 ) have 
the same value. 

demo file for target option FROM_FILE (homogeneous, isotropic target) 



111) 

of "target origin' 



8 








NAT 








1 . 


000 




0, 


. 000 


. 


000 


= target vector al (in TF) 


0. 


000 




1 


. 000 


. 


000 


= target vector a2 (in TF) 


1 . 






1 




1 . 




= d_x/d d_y/d d_z/d (normally 


0. 


5 




0. 


,5 


0. 


5 


= X0(l-3) = location in lattice 


J 




JX 




JY 


JZ 


ICOMPX, ICOMPY, ICOMPZ 


1 















111 


2 












1 


111 


3 









1 





111 


4 









1 


1 


111 


5 




1 
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6 




1 







1 


111 


7 




1 




1 





111 


8 




1 




1 


1 


111 



The above sample target consists of 8 dipoles arranged to represent a cube. 
This example is homogeneous: All sites have composition 1 
The target origin X is set to be at the center of the target 

Note that ICOMPX,ICOMPY,ICOMPZ could differ, allowing treatment of anisotropic targets, provided 
the dielectric tensor at each location is diagonal in the TF. 

demo file for target option FROM_FILE ( inhomogeneous , isotropic target) — 
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8 








NAT 








1 . 


000 




0. 


. 000 


0. 


000 


= target vector al (in TF) 


0. 


000 




1. 


. 000 


0. 


000 


= target vector a2 (in TF) 


1 . 






1 , 




1 . 




= d_x/d d_y/d d_z/d (normally 


0. 


5 




0, 


.5 


0. 


5 


= X0(l-3) = location in lattice 


J 




JX 




JY 


JZ 


ICOMPX, ICOMPY, ICOMPZ 


1 
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1 





111 


4 









1 


1 


111 


5 




1 










2 2 2 


6 




1 







1 


2 2 2 


7 




1 




1 





2 2 2 


8 




1 




1 


1 


2 2 2 



TARGET GENERATION: ISOLATED FINITE TARGETS 



This sample target consists of 8 dipoles arranged to represent a cube. 

This example is inhomogeneous: The lower half of the cube (JX=0) has isotropic composition 1 
The upper half of the cube (JX=1) has isotropic composition 2 
The target origin XO is set to be at the center of the target. 

Note that ICOMPX, ICOMPY, ICOMPZ can be different, allowing treatment of anisotropic targets, 
provided the dielectric tensor at each location is diagonal in the TF. 

Sample calculation: The directory examples_exp/FROM_FILE contains files for calculation of 
scattering by a 0.5/im x 1/im x Ifim Au block, represented by a 32 x 64 x 64 dipole array. 



19.2 ANIFRMFIL = General anistropic target denned by list of dipole loca- 
tions,"compositions", and material orientations obtained from a file 

This option causes DDSCAT to read the target geometry information from a file shape . dat instead 
of automatically generating one of the geometries listed below. The file shape . dat gives the num- 
ber N of dipoles in the target, the components of the "target axes" ai and a2 in the Target Frame 
(TF), the vector xq(1 — 3) determining the correspondence between the integers IXYZ and actual co- 
ordinates in the TF, and specifications for the location and "composition" of each dipole. For each 
dipole J, the file shape . dat provides the location IXYZ ( J, 1-3 ) , the composition identifier integer 
I COMP (J, 1 - 3 ) specifying the dielectric function corresponding to the three principal axes of the di- 
electric tensor, and angles 8df> $df. an d /?df specifying the orientation of the local "Dielectric Frame" 
(DF) relative to the "Target Frame" (TF) (see $26]). The DF is the reference frame in which the dielectric 
tensor is diagonalized. The Target Frame is the reference frame in which we specify the dipole locations. 

The shape . dat file is read by routine REASHP (file reashp . f 90). The user can customize 
REASHP as needed to conform to the manner in which the target geometry is stored in file shape . dat. 
However, as supplied, REASHP expects the file shape . dat to have the following structure: 

• one line containing a description; the first 67 characters will be read and printed in various output 
statements. 

• N = number of dipoles in target 

• CLix &iy a iz = x,y,z components (in Target Frame) of ai 

• <J,2x ci2y a-2z = x,y,z components (in Target Frame) of &2 

• d x /d dy/d d z /d = 1. 1. 1. = relative spacing of dipoles in xtf, Ytf, ztf directions 

• xox xo y xqz = TF coordinates xtf / d yrp/d zrp/d corresponding to lattice site IXYZ= 

• (line containing comments) 

• dummy IXYZ (1, 1-3) ICOMP (1,1-3) THETADF ( 1 ) PHIDF(l) BETADF ( 1 ) 

• dummy IXYZ (2, 1-3) ICOMP (2, 1-3) THETADF (2) PHIDF(2) BETADF (2) 

• dummy IXYZ (3, 1-3) ICOMP (3, 1-3) THETADF (3) PHIDF(3) BETADF (3) 



dummy IXYZ (J , 1-3) ICOMP (J, 1-3) THETADF (J) PHIDF(J) BETADF (J) 



• dummy IXYZ (N, 1-3) ICOMP (N, 1-3) THETADF (N) PHIDF (N) BETADF (N) 

Where dummy is a number (either integer or floating point) that might, for example, give the number 
identifying the dipole. This number will not be used in any calculations. 
THETADF PHIDF BETADF should be given in radians. 

Here is an example of the first few lines of a target description file suitable for target option AN I F RMF I L 

demo file for target option ANIFRMFIL (this line is for comments) 

8 = NAT 

1.000 0.000 0.000 = target vector al (in TF) 
0.000 1.000 0.000 = target vector a2 (in TF) 
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This sample target consists of 8 dipoles arranged to represent a cube. 
Half of the cube (dipoles with JX=0) has isotropic composition 1 . For this case, the angles THETADF, 
PHIDF, BETADF do not matter, and it convenient to set them all to zero. 

The other half of the cube (dipoles with JX=1) consists of a uniaxial material, with dielectric function 2 
for E fields parallel to one axis (the "c-axis"), and dielectric function 3 for E fields perpendicular to the 
c-axis. The c-axis is 30° (0.5236 radians) away from xtf, and lies in the xtf-ztf plane (having been 
rotated by 1.5708 radians around Xtf)- 

Note that ICOMP(J,K) can be different for K=l,3, allowing treatment of anisotropic targets, provided 
the dielectric tensor at each location is diagonal in the TF. 

19.3 ANIELLIPS = Homogeneous, anisotropic ellipsoid. 

SHPARi, SHPAR 2 , SHPAR 3 define the ellipsoidal boundary: 

/ XTF/d \ 2 ( y TF /d \ 2 / z TF /d \ 2 = 1 
V SHPARi/ V SHPAR 2/ \SHPAR 3 / 4 ' 

The TF origin is located at the centroid of the ellipsoid. 



19.4 ANI_ELL_2 = Two touching, homogeneous, anisotropic ellipsoids, with dis- 
tinct compositions 

Geometry as for ELLIPS0.2; SHPARi, SHPAR 2 , SHPAR 3 have same meanings as for ELLIPS0.2. 

Target axes ai = (1, 0, 0)tf and a 2 = (0,1, 0)tf- 

Line connecting ellipsoid centers is || ai = xtf- 

TF origin is located between ellipsoids, at point of contact. 

It is assumed that (for both ellipsoids) the dielectric tensor is diagonal in the TF. User must set NCOMP=6 
and provide xx, yy, zz components of dielectric tensor for first ellipsoid, and xx, yy, zz components of 
dielectric tensor for second ellipsoid (ellipsoids are in order of increasing xtf)- 



19.5 ANI_ELL_3 = Three touching homogeneous, anisotropic ellipsoids with 
same size and orientation but distinct dielectric tensors 

SHPARi, SHPAR 2 , SHPAR 3 have same meanings as for ELLIPS0.3. 

Target axis ai = (1, 0, 0)tf (along line of ellipsoid centers), and a 2 = (0, 1, 0)tf- 

TF origin is located at center of middle ellipsoid. 

It is assumed that dielectric tensors are all diagonal in the TF. User must set NCOMP=9 and provide xx, 
yy, zz elements of dielectric tensor for first ellipsoid, xx, yy, zz elements for second ellipsoid, and xx, 
yy, zz elements for third ellipsoid (ellipsoids are in order of increasing xtf)- 



19.6 ANIRCTNGL = Homogeneous, anisotropic, rectangular solid 

x, y, z lengths/rf= SHPARi, SHPAR 2 , SHPAR 3 . 

Target axes ai = (1,0, 0)tf and a 2 = (0,1, 0)tf in the TF. 
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(xtf, Utf, ztf) = (0, 0, 0) at middle of upper target surface, (where "up" = xtf)- (The target surface 

is taken to be d/2 about the upper dipole layer.) 

Dielectric tensor is assumed to be diagonal in the target frame. 

User must set NCOMP=3 and supply names of three files for e as a function of wavelength or energy: first 
for e xx , second for e yy , and third for e zz , as in the following sample ddscat . par file: 

' ============ Parameter file for v7 . 1 . ' 

'**** PRELIMINARIES 

'NOTORQ'= CMTORQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 

'PBCGS2'= CMDSOL*6 (PBCGS2, PBCGST, PETRKP) — select solution method 

' GPFAFT' = CMETHD*6 (GPFAFT, FFTMKL) 

' GKDLDR' = CALPHA*6 (GKDLDR, LATTDR) 

'NOTBIN'= CBINFLAG (ALLBIN, ORIBIN, NOTBIN) 

' ANIRCTNGL' = CSHAPE*9 shape directive 

'**** Initial Memory Allocation ****' 

10 25 50 = upper bound on target extent 

'**** Target Geometry and Composition 

' ANIRCTNGL' = CSHAPE*9 shape directive 

10 25 50 = shape parameters SHPAR1, SHPAR2, SHPAR3 
3 = NCOMP = number of dielectric materials 

' /u/draine/DDA/diel/m2 . 0_0 . 1 ' = name of file containing dielectric function 
' /u/draine/DDA/diel/ml . 5 0„0 .00' 
' /u/draine/DDA/diel/ml . 5 0„0 .00' 
'**** Error Tolerance ****' 

1.0 0e-5 = TOL = MAX ALLOWED (NORM OF | G>=AC I E>-ACA | X> ) / (NORM OF AC|E>) 

'**** Interaction cutoff parameter for PBC calculations ****' 

5.00e-3 = GAMMA (le-2 is normal, 3e-3 for greater accuracy) 

'**** Angular resolution for calculation of <cos>, etc. ****' 

2.0 = ETASCA (number of angles is proportional to [ (2+x) /ETASCA] "2 ) 

'**** Wavelengths (micron) ****' 

500. 500. 1 ' INV = wavelengths ( first , last , how many , how=LIN, INV, LOG) 
' **** Effective Radii (micron) **** ' 

30.60 30.60 1 'LIN' = eff. radii (first, last, how many, how=LIN, INV, LOG) 
'**** Define Incident Polarizations ****' 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = IORTH (=1 to do only pol . state eOl; =2 to also do orth. pol . state) 

'**** Specify which output files to write ****' 

1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient. 

1 = IWRPOL (=0 to suppress, =1 to write ".pol" file for each (BETA, THETA) 
' **** Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 
0. 90. 3 = THETMI, THETMX, NTHETA (theta=angle between al and k) 
0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 
'**** Specify first IWAV, IRAD, IORI (normally 0) ****' 
= first IWAV, first IRAD, first IORI (0 to begin fresh) 
' **** Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 

11 12 21 22 31 41 = indices ij of elements to print 
'**** Specify Scattered Directions ****' 

' LFRAME' = CMDFRM* 6 ( ' LFRAME ' or ' TFRAME ' for Lab Frame or Target Frame) 

2 = number of scattering planes 

0. 0. 180. 30 = phi, thetan_min, thetan_max, dtheta (in degrees) for plane A 
90. 0. 180. 30 = phi, ... for plane B 



19.7 CONELLIPS = Two concentric ellipsoids 

SHPARi, SHPAR 2 , SHPAR 3 = lengths/d of the outer ellipsoid along the x T f, Ytf, z T f axes; 
SHPAR4, SHPAR5, SHPAR6 = lengths/c? of the inner ellipsoid along the xtf, Ytf, ztf axes. 
Target axes ai = (1, 0, 0)tf, a 2 = (0, 1, 0)tf- 
TF origin is located at centroids of ellipsoids. 

The "core" within the inner ellipsoid is composed of isotropic material 1 ; the "mantle" between inner 
and outer ellipsoids is composed of isotropic material 2. 

User must set NCOMP=2 and provide dielectric functions for "core" and "mantle" materials. 
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19.8 CYLINDER1 = Homogeneous, isotropic finite cylinder 

SHPARi = length/d, SHPAR 2 = diameter/d, with 

SHPAR3 = 1 for cylinder axis ai || xtf: ai = (1, 0, 0)tf and a 2 = (0, 1, 0)tf; 
SHPAR3 = 2 for cylinder axis ai || Ytf: &i = (0,1, 0)tf and a 2 = (0, 0, 1)tf; 
SHPAR3 = 3 for cylinder axis ai || ztf: &i = (0, 0, 1)tf and a 2 = (1, 0, 0)tf in the TF. 
TF origin is located at centroid of cylinder. 
User must set NCOMP=l. 



19.9 CYLNDRCAP = Homogeneous, isotropic finite cylinder with hemispheri- 
cal endcaps. 

SHPARi = cylinder length/d (not including end-caps!) and SHPAR 2 = cylinder diameter/d, with cylinder 

axis = ai = (1,0, 0)tf and a 2 = (0, 1,0)tf- The total length along the target axis (including the 

endcaps) is (SHPARi+SHPAR 2 )d. 

TF origin is located at centroid of cylinder. 

User must set NCOMP=l. 



19.10 DSKRCTNGL = Disk on top of a homogeneous rectangular slab 

This option causes DDSCAT to create a target consisting of a disk of composition 1 resting on top of a 
rectangular block of composition 2. Materials 1 and 2 are assumed to be homogeneous and isotropic, 
ddscat . par should set NCOMP to 2 . 

The cylindrical disk has thickness SHPARi xd in the x-direction, and diameter SHPAR 2 xd. The rect- 
angular block is assumed to have thickness SHPARaxd in the x-direction, length SHPAR4xd in the 
y-direction, and length SHPARsxd in the z-direction. The lower surface of the cylindrical disk is in the 
x = plane. The upper surface of the slab is also in the x = plane. 

The Target Frame origin (0,0,0) is located where the symmetry axis of the disk intersects the x = 
plane (the upper surface of the slab, and the lower surface of the disk). In the Target Frame, dipoles 
representing the rectangular block are located at (x/d, y/d, z/d) = (j x + 0.5, j y + A y ,j z + A z ), where 
jx, jy, and j z are integers. A y = or 0.5 depending on whether SHPAR4 is even or odd. A 2 = or 0.5 
depending on whether SHPAR5 is even or odd. 

Dipoles representing the disk are located at 
x/d = 0.5, 1.5, [int(SHPAR 4 + 0.5) - 0.5] 

As always, the physical size of the target is fixed by specifying the value of the effective radius 
a c ff = (3VT/47T) 1 / 3 , where Vr is the total volume of solid material in the target. For this geome- 
try, the number of dipoles in the target will be approximately N = [SHPARi x SHPAR 2 x SHPAR3 + 
(7r/4)((SHPAR4) 2 x SHPAR5)], although the exact number may differ because of the dipoles are required 
to be located on a rectangular lattice. The dipole spacing d in physical units is determined from the 
specified value of a e ff and the number N of dipoles in the target: d = (47r/3iV) 1 / 3 a ff . This option 
requires 5 shape parameters: 
The pertinent line in ddscat . par should read 



SHPARI SHPAR 2 SHPAR3 SHPAR4 SHPAR5 



where 

SHPARi = [disk thickness (in xtf direction)]/^ [material 1] 
SHPAR 2 = (disk diameter)/d 

SHPAR3 = (brick thickness in xtf direction)/^ [material 2] 
SHPAR4 = (brick thickness in yxF direction)/d 
SHPAR5 = (brick thickness in ztf direction)/d 



The overall size of the target (in terms of numbers of dipoles) is determined by parameters ( SHPARI + SHPAR4 ) , 
SHPAR 2 , and SHPAR3. The periodicity in the TF y and z directions is determined by parameters SHPAR4 
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and SHPAR5. 

The physical size of the TUC is specified by the value of a e fi (in physical units, e.g. cm), specified in the 
file ddscat . par. 

The "computational volume" is determined by (SHPARi + SHPAR 4 )x SHPAR 2 x SHPAR 3 . 

The target axes (in the TF) are set to ai = xtf = (1, 0, 0)tf - i-e., normal to the "slab" - and 
&2 = Ytf = (0,1, 0)tf- The orientation of the incident radiation relative to the target is, as for all other 
targets, set by the usual orientation angles 9, $, and j3 (see SjTT] above); for example, O = would be 
for radiation incident normal to the slab. 

19.11 DW1996TAR = 13 block target used by Draine & Weingartner (1996). 

Single, isotropic material. Target geometry was used in study by Draine & Weingartner (1996) of radia- 
tive torques on irregular grains, ai and a 2 are principal axes with largest and second-largest moments of 
inertia. User must set NCOMP=l. Target size is controlled by shape parameter SHPAR ( 1 ) = width of 
one block in lattice units. 
TF origin is located at centroid of target. 

19.12 ELLIPSOID = Homogeneous, isotropic ellipsoid. 

"Lengths" SHPARi, SHPAR 2 , SHPAR 3 in the x, y, z directions in the TF: 

( - ttf V + ( V+ ( ztf V = I ( 33) 

VsHPAMy V SHPAR 2<v V SHPAR 3 rf / 4 ' 

where d is the interdipole spacing. 

The target axes are set to ai = (1, 0, 0)tf and a 2 = (0, 1, 0)tf- 

Target Frame origin = centroid of ellipsoid. 

User must set NCOMP=l on line 9 of ddscat . par. 

A homogeneous, isotropic sphere is obtained by setting SHPARi = SHPAR 2 = SHPAR3 = diameter/rf. 

19.13 ELLIPSO_2 = Two touching, homogeneous, isotropic ellipsoids, with dis- 
tinct compositions 

SHPARi, SHPAR 2 , SHPAR 3 =x-length/d, y-length/<i, z-lmgth/d of one ellipsoid. The two ellipsoids 
have identical shape, size, and orientation, but distinct dielectric functions. The line connecting ellipsoid 
centers is along the XTF-axis. Target axes ai = (1,0,0)tf [along line connecting ellipsoids] and 
a 2 = (0,1,0) TF . 

Target Frame origin = midpoint between ellipsoids (where ellipsoids touch). 

User must set NCOMP=2 and provide dielectric function file names for both ellipsoids. Ellipsoids are 
in order of increasing x^p: first dielectric function is for ellipsoid with center at negative xtf, second 
dielectric function for ellipsoid with center at positive xtf • 

19.14 ELLIPSO_3 = Three touching homogeneous, isotropic ellipsoids of equal 
size and orientation, but distinct compositions 

SHPARi, SHPAR 2 , SHPAR3 have same meaning as forELLIPS0_2. Line connecting ellipsoid centers is 
parallel to xtf axis. Target axis ai = (1, 0, 0)tf (along line of ellipsoid centers), and a 2 = (0, 1, 0)tf- 
Target Frame origin = centroid of middle ellipsoid. 

User must set NCOMP=3 and provide (isotropic) dielectric functions for first, second, and third ellipsoid. 

Sample calculation: The directory examples_exp/ELLIPSOID contains ddscat . par for cal- 
culation of scattering by a sphere with refractive index m = 1.5 + O.Oli and 27ra/A = 5, represented by 
a N = 110592 dipole pseudosphere just fitting within a 48 x 48 x 48 computational volume. 
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19.15 HEX_PRISM = Homogeneous, isotropic hexagonal prism 

SHPARi = (Length of prism)/d = (distance between hexagonal faces)/d, 

SHPAR2 = (distance between opposite vertices of one hexagonal face)/e? = 2 x hexagon side/eL 

SHPAR3 selects one of 6 orientations of the prism in the Target Frame (TF). 

Target axis a.i is along the prism axis (i.e., normal to the hexagonal faces), and target axis a 2 is normal 
to one of the rectangular faces. There are 6 options for SHPAR3: 

SHPAR3 = 1 for a! || Xtf and a 2 || ytf ; SHPAR3 = 2 for a.! || xtf an d a 2 II ^tf ; 
SHPAR3 = 3 for a! || ytf an d &2 II ^tf ; SHPAR3 = 4 for a.! || yxF an d a 2 II ^tf ; 
SHPAR3 = 5 for a! || Zxf an d a 2 || Xtf ; SHPAR 3 = 6 for a.! || Ztf an d a 2 || Ytf 
TF origin is located at the centroid of the target. 
User must set NC0MP=1. 

19.16 LAYRDSLAB = Multilayer rectangular slab 

Multilayer rectangular slab with overall x, y, z lengths a x = SHPARi x d 
a y = SHPAR2 x d, 
a z = SHPAR3 x d. 

Upper surface is at xtf = 0, lower surface at xtf = —SHPARi x d 

SHPAR4 = fraction which is composition 1 (top layer). 

SHPAR5 = fraction which is composition 2 (layer below top) 

SHPAR6 = fraction which is composition 3 (layer below comp 2) 

1 - (SHPAR4 + SHPAR 5 + SHPAR 6 ) = fraction which is composition 4 (bottom layer). 

To create a bilayer slab, just set SHPAR 5 = SHPAR 6 = 

To create a trilayer slab, just set SHPAR6 = 

User must set NCOMP=2,3, or 4 and provide dielectric function files for each of the two layers. Top 
dipole layer is at xtf = —d/2. Origin of TF is at center of top surface. 

19.17 MLTBLOCKS = Homogeneous target constructed from cubic "blocks" 

Number and location of blocks are specified in separate file blocks. par with following structure: 
one line of comments (may be blank) 
PRIN (= or 1 - see below) 
N (= number of blocks) 
B (= width/tf of one block) 

xtf Utf ztf (= position of 1st block in units of Bd) 
Xtf Utf ztf (= position of 2nd block in units of Bd) ) 

xtf Utf ztf (= position of Nth block in units of Bd) 
If PRIN=0, then ai = (1, 0, 0)tf, a 2 = (0, 1, 0)tf- If PRIN=1, then ai and a 2 are set to principal axes 
with largest and second largest moments of inertia, assuming target to be of uniform density. User must 

set NCOMP=l. 

19.18 RCTGLPRSM = Homogeneous, isotropic, rectangular solid 

x, y, z lengths/d = SHPARi, SHPAR 2 , SHPAR3. 
Target axes ai = (1,0, 0)tf an d a 2 = (0, 1, 0)tf- 

TF origin at center of upper surface of solid: target extends from xtf /d = —SHPARi to 0, 
y T F/d from -0.5 x SHPAR 2 to +0.5 x SHPAR 2 
ZTF/d from -0.5 x SHPAR 3 to +0.5 x SHPAR 3 
User must set NCOMP=l. 

The directory examples_exp/RCTGLPRSM contains ddscat .par for calculation of scattering 
by a 0.5/^m x l[im x lfim Au block, represented by a 32 x 64 x 64 dipole array. 
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19.19 

axis. 



RCTGLBLK3 = Stack of 3 rectangular blocks, with centers on the x TF 



Each block consists of a distinct material. There are 9 shape parameters: 

SHPARi = (upper solid thickness in xtf direction)/d [material 1] 

SHPAR2 = (upper solid width in yxF direction)/^ 

SHPAR3 = (upper solid width in ztf direction)/^ 

SHPAR4 = (middle solid thickness in xtf direction)/d [material 2] 

SHPAR5 = (middle solid width in yxF direction)/d 

SHPARg = (middle solid width in ztf direction)/d 

SHPAR7 = (lower solid thickness in xtf direction)/d [material 3] 

SHPARg = (lower solid width in yxF direction)/^ 

SHPARg = (lower solid width in ztf direction)/^ 

TF origin is at center of top surface of material 1 . 



19.20 SLAB_HOLE = Rectangular slab with a cylindrical hole. 

The target consists of a rectangular block with a cylindrical hole with the axis passing through the 
centroid and aligned with the xtf axis. The block dimensions are a x b x c. The cylindrical hole has 
radius r. The pertinent line in ddscat . par should read 

SHPARi SHPAR 2 SHPAR3 SHPAR4 
where 

SHPARi = a/d(d is the interdipole spacing) 

SHPAR 2 = b/a 
SHPAR3 = c/a 
SHPAR4 = r/a 

Ideally, SHPARi, SHPAR 2 x SHPARi, SHPAR 3 x SHPARi will be integers (so that the cubic lattice can 
accurately approximate the desired target geometry), and SHPAR4X SHPARi will be large enough for 
the circular cross section to be well-approximated. 

The TF origin is at the center of the top surface (the top surface lies in the yxF — ztf plane, and 
extends from i/tf = —b/2 to +6/2, and ztf = — c/2 to +c/2). The cylindrical hole axis runs from 
{xtf = 0, ?/tf = 0, z TF = 0) to (x T f = -a, 2/tf = 0, ztf = 0). 



19.21 SPHERES_N = Multisphere target = union of N spheres of single isotropic 
material 

Spheres may overlap if desired. The relative locations and sizes of these spheres are defined in an exter- 
nal file, whose name (enclosed in single quotes) is passed through ddscat . par. The length of the file 
name should not exceed 80 characters. The pertinent line in ddscat . par should read 

SHPARi SHPAR2 'filename ' (quotes must be used) 

where SHPARi = target diameter in x direction (in Target Frame) in units of d 
SHPAR 2 = to have a\ = (1, 0, 0) T f, a 2 = (0, 1, 0) T f- 

SHPAR2= 1 to use principal axes of moment of inertia tensor for a± (largest /) and a 2 (intermediate /). 
filename is the name of the file specifying the locations and relative sizes of the spheres. 
The overall size of the multisphere target (in terms of numbers of dipoles) is determined by parameter 
SHPARi, which is the extent of the multisphere target in the x-direction, in units of the lattice spacing d. 
The file filename ' should have the following structure: 



N (= number of spheres) 

line of comments (may be blank) 

line of comments (may be blank) [N.B.: changed from v7.0.7] 
line of comments (may be blank) [N.B.: changed from v7.0.7] 
line of comments (may be blank) [N.B.: changed from v7.0.7] 
x\ yi z\ a\ (arb. units) 
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x 2 2/2 z 2 a 2 (arb. units) 
xn Un zn o>n (arb. units) 

where Xj, yj, Zj are the coordinates (in the TF) of the center of sphere j, and aj is the radius of sphere j. 
Note that Xj, yj, Zj, dj (j — 1, N) establish only the shape of the N— sphere target. For instance, 
a target consisting of two touching spheres with the line between centers parallel to the x axis could 
equally well be described by lines 6 and 7 being 

0.5 

1 0.5 

or 

1 

2 1 

The actual size (in physical units) is set by the value of a e ff specified in ddscat . par, where, as always, 
a c ff = (3V/47T) 1 / 3 , where V is the total volume of material in the target. 

Target axes ai and a 2 are set to be principal axes of moment of inertia tensor (for uniform density), 
where ai corresponds to the largest eigenvalue, and a. 2 to the intermediate eigenvalue. 
The TF origin is taken to be located at the volume-weighted centroid. 
User must set NCOMP=l. 

19.22 SPHROID_2 = Two touching homogeneous, isotropic spheroids, with dis- 
tinct compositions 

First spheroid has length SHPARi along symmetry axis, diameter SHPAR2 perpendicular to symmetry 
axis. Second spheroid has length SHPAR3 along symmetry axis, diameter SHPAR4 perpendicular to 
symmetry axis. Contact point is on line connecting centroids. Line connecting centroids is in xtf 
direction. Symmetry axis of first spheroid is in yTF direction. Symmetry axis of second spheroid is in 
direction y T F cos(SHPAR 5 ) + z T f sin(SHPAR 5 ), and SHPAR5 is in degrees. If SHPAR 6 = 0., then target 
axes ai = (1, 0, 0)tf> &2 = (0, 1, 0)tf- If SHPARg = 1., then axes ai and a2 are set to principal axes 
with largest and 2nd largest moments of inertia assuming spheroids to be of uniform density. 
Origin of TF is located between spheroids, at point of contact. 
User must set NCOMP=2 and provide dielectric function files for each spheroid. 

19.23 SPH_ANI_N = Multisphere target consisting of the union of TV spheres of 
various materials, possibly anisotropic 

Spheres may NOT overlap. The relative locations and sizes of these spheres are defined in an external 
file, whose name (enclosed in single quotes) is passed through ddscat .par. The length of the file 
name should not exceed 80 characters. Target axes ai and & 2 are set to be principal axes of moment 
of inertia tensor (for uniform density), where ai corresponds to the largest eigenvalue, and a.2 to the 
intermediate eigenvalue. 

The TF origin is taken to be located at the volume-weighted centroid. 
The pertinent line in ddscat . par should read 

SHPARi SHPAR2 'filename ' (quotes must be used) 

where SHPARi = target diameter in x direction (in Target Frame) in units of d 
SHPAR 2 = to have ai = (1, 0, 0)tf, a 2 = (0, 1, 0) T f in Target Frame. 

SHPAR2= 1 to use principal axes of moment of inertia tensor for a\ (largest /) and a 2 (intermediate /). 
filename is the name of the file specifying the locations and relative sizes of the spheres. 
The overall size of the multisphere target (in terms of numbers of dipoles) is determined by parameter 
SHPARi, which is the extent of the multisphere target in the x-direction, in units of the lattice spacing d. 
The file filename ' should have the following structure: 

TV (= number of spheres) 

line of comments (may be blank) 
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line of comments (may be blank) [N.B.: changed from v7.0.7] 
line of comments (may be blank) [N.B.: changed from v7.0.7] 
line of comments (may be blank) [N.B.: changed from v7.0.7] 

xi yi zi n Cxi Cyi Cz\ 0df,i $df,i /Sdf.i 
x 2 2/2 z 2 r 2 Cx 2 Cy 2 Cz 2 6 D f,2 $df,2 /?df,2 

xn Vn z n r N Cx N Cy N Cz N Qbf.n $df,n Pdf,n 

where Xj, jjj , Zj £irc the coordinates of the center, and rj is the radius of sphere j (arbitrary units), 
Cxj, Cyj, Czj are integers specifying the "composition" of sphere j in the x,y,z directions in the 
"Dielectric Frame" (see £|26l i of sphere j, and ©dfj ^ dfj Pdfj are angles (in radians) specifying 
orientation of the dielectric frame (DF) of sphere j relative to the Target Frame. Note that Xj, yj, zj, rj 
(j = 1, N) establish only the shape of the N— sphere target, just as for target option NSPHER. The 
actual size (in physical units) is set by the value of a g specified in ddscat . par, where, as always, 
a e ff = (SV/Att) 1 / 3 , where V is the volume of material in the target. 

User must set NCOMP to the number of different dielectric functions being invoked (i.e., the range of 

{Cxj,Cyj,Czj}. 

Note that while the spheres can be anisotropic and of differing composition, they can of course also 
be isotropic and of a single composition, in which case the relevant lines in the file 'filename ' would be 
simply 

iV (= number of spheres) 

line of comments (may be blank) 

line of comments (may be blank) [N.B.: changed from v7.0.7] 
line of comments (may be blank) [N.B.: changed from v7.0.7] 
line of comments (may be blank) [N.B.: changed from v7.0.7] 
X! yi zi n 1 1 1 
x 2 V2 z 2 r 2 1 1 1 

xn Un zn rN 1 1 1 

Subdirectory examples_exp/SPH_ANI_N contains ddscat .par for calculating scattering by 
a random aggregate of 64 spheres [aggregated according following the "BAM2" aggregation process 
described by Shen, Draine & Johnson (2008)]. 32 of the spheres are assumed to consist of "astrosilicate", 
and 32 of crystalline graphite, with random orientations for each of the 32 graphite spheres. Each sphere 
is assumed to have a radius 0.050^m. The entire cluster is represented by N = 7947 dipoles, or about 
124 dipoles per sphere. Scattering and absorption are calculated for A = 0.55/im. 

19.24 TETRAHDRN = Homogeneous, isotropic tetrahedron 

SHPARi=length/<i of one edge. Orientation: one face parallel to yTF, ztf plane, opposite "vertex" is 
in +xtf direction, and one edge is parallel to ztf- Target axes ai = (1, 0, 0)tf [emerging from one 
vertex] and & 2 — (0, 1, 0)tf [emerging from an edge] in the TF. User must set NCOMP=l. 

19.25 TRNGLPRSM = Triangular prism of homogeneous, isotropic material 

SHPARi, SHPAR 2 , SHPAR 3 , SHPAR 4 = a/d, b/a, c/a, L/a 

The triangular cross section has sides of width a, b, c. L is the length of the prism, d is the lattice 
spacing. The triangular cross-section has interior angles a, (3, 7 (opposite sides a, b, c) given by cos a = 
(b 2 + c 2 - a 2 )/2bc, cos/3 = (a 2 + c 2 - b 2 )/2ac, cos 7 = (a 2 + b 2 - c 2 )/2ab. In the Target Frame, 
the prism axis is in the x direction, the normal to the rectangular face of width a is (0,1,0), the normal to 
the rectangular face of width b is (0, — cos 7, sin 7), and the normal to the rectangular face of width c is 
(0, — cos 7, — sin 7). 
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19.26 UNIAXICYL = Homogeneous finite cylinder with uniaxial anisotropic di- 
electric tensor 

SHPARi, SHPAR 2 have same meaning as for CYLINDER].. Cylinder axis = ai = (1,0, 0)tf, a 2 = 
(0, 1, 0)tf- It is assumed that the dielectric tensor e is diagonal in the TF, with e yy = e zz . User must set 
NCOMP=2. Dielectric function 1 is for E || ai (cylinder axis), dielectric function 2 is for E _L a.\. 



19.27 Modifying Existing Routines or Writing New Ones 

The user should be able to easily modify these routines, or write new routines, to generate targets with 
other geometries. The user should first examine the routine target . f 90 and modify it to call any 
new target generation routines desired. Alternatively, targets may be generated separately, and the target 
description (locations of dipoles and "composition" corresponding to x,y,z dielectric properties at each 
dipole site) read in from a file by invoking the option FROM_FILE in ddscat . f 90. 

Note that it will also be necessary to modify the routine reapar . f 90 so that it will accept whatever 
new target option is added to the target generation code . 



19.28 Testing Target Generation using CALLTARGET 

It is often desirable to be able to run the target generation routines without running the entire DDSCAT 
code. We have therefore provided a program CALLTARGET which allows the user to generate targets 
interactively; to create this executable just type^l 

make calltarget . 

The program calltarget is to be run interactively; the prompts are self-explanatory. You may need 
to edit the code to change the device number IDVOUT as for DDSCAT (see %J\ abo\e). 

After running, calltarget will leave behind an ASCII file target . out which is a list of the 
occupied lattice sites in the last target generated. The format of target . out is the same as the format 
of the shape . dat files read if option FROM_FILE is used (see above). Therefore you can simply 

mv target. out shape.dat 
and then use DDSCAT with the option FROM_FILE (or option ANIFRMFIL in the case of anisotropic 
target materials with arbitrary orientation relative to the Target Frame) in order to input a target shape 
generated by CALLTARGET. 

Note that CALLTARGET - designed to generate finite targets - can be used with some of the "PBC" 
target options (see i}20l below) to generate a list of dipoles in the TUC. At the moment, CALLTARGET 
has support for target options BISLINPBC, DSKBLYPBC, and DSKRCTPBC. 

20 Target Generation: Periodic Targets 

A periodic target consists of a "Target Unit Cell" (TUC) which is then repeated in either the yTF direc- 
tion, the ztf direction, or both. Please see Draine & Flatau (2008) for illustration of how periodic targets 
are assembled out of TUCs, and how the scattering from these targets in different diffraction orders M 
or (M, N) is constrained by the periodicity. 

The following options are included in DDSCAT: 

9 Non-Linux sites: The source code for CALLTARGET is in the file CALLTARGET . f 90. You must compile and link 

CALLTARGET . f 90, dclcommon . f 90, dsyevj3.f90, errmsg.f90, gasdev.f90, p.lm.f90, prinaxis . f 90, 
ran3.f90, reashp.f90, sizer.f90, tar2el.f90, tar2sp.f90, tar3el.f90, taranirec . f 90, 
tarblocks . f 90, tarcel.f90, tarcyl.f90, tarcylcap . f 90, tarell.f90, target. f90, targspher . f 90, 
tarhex.f90, tarnas.f90, tarnsp.f90, tarpbxn.f90, tarprsm. f 90, tarrctblk3 . f 90, tarrecrec . f 90, 
tarslblin . f 90, tartet . f 90, and wrimsg . f 90. 
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20.1 FRMFILPBC = periodic target with TUC geometry and composition input 
from a file 

The TUC can have arbitrary geometry and inhomogeneous composition, and is assumed to repeat peri- 
odically in either 1-d (y or z) or 2-d (y and z). 
The pertinent line in dd scat. par should read 
SHPARi SHPAR2 'filename ' (quotes must be used) 
SHPARi = Py/d (P y = periodicity in yxF direction) 
SHPAR2 = P z /d (P z = periodicity in ztf direction) 

filename is the name of the file specifying the locations of the dipoles, and the "composition" at each 
dipole location. The composition can be anisotropic, but the dielectric tensor must be diagonal in the 
TF. The shape and composition of the TUC are provided exactly as for target option FROM_FILE - see 

em 

If SHPARi = then the target does not repeat in the yxF direction. 
If SHPAR2 = then the target does not repeat in the ztf direction. 

20.2 ANIFILPBC = general anisotropic periodic target with TUC geometry and 
composition input from a file 

The TUC can have arbitrary geometry and inhomogeneous composition, and is assumed to repeat peri- 
odically in either 1-d (y or z) or 2-d (y and z). 
The pertinent line in dd scat. par should read 
SHPARi SHPAR2 'filename ' (quotes must be used) 
SHPARi = Py/d (P y = periodicity in yxF direction) 
SHPAR2 = P z /d (P z = periodicity in ztf direction) 

filename is the name of the file specifying the locations of the dipoles, and the "composition" at each 
dipole location. The composition can be anisotropic, and the dielectric tensor need not be diagonal in 
the TF. The shape and composition of the TUC are provided exactly as for target option ANIFRMF IL - 
see COO] 

If SHPARi = then the target does not repeat in the yxF direction. 
If SHPAR2 = then the target does not repeat in the ztf direction. 



20.3 BISLINPBC = Bi-Layer Slab with Parallel Lines 

The target consists of a bi-layer slab, on top of which there is a "line" with rectangular cross-section. 
The "line" on top is composed of material 1, has height X\ (in the xtf direction), width Y\ (in the yxF 
direction), and is infinite in extent in the ztf direction. 

The bilayer slab has width Y2 (in the yxF direction). It is consists of a layer of thickness X 2 of material 
2, on top of a layer of material 3 with thickness X3. 



SHPARi 


= X x /d 


(Xi 


= thickness of line) 


SHPAR 2 


= Y 1 /d 


(Yi = 


: width of line) 


SHPAR3 


= X 2 /d 


(X 2 


= thickness of upper layer of slab) 


SHPAR4 


= X 3 /d 


(X 3 


= thickness of lower layer of slab) 


SHPAR5 


= Y 2 /d 


(Y 2 = 


: width of slab) 


SHPAR 6 


= Py/d 


(Py- 


= periodicity in yxF direction). 



If SHPARg = 0, the target is NOT periodic in the yxF direction, consisting of a single column, infinite 
in the ztf direction. 
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20.4 CYLNDRPBC = Target consisting of homogeneous cylinder repeated in 
target y and/or z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of an infinite array of cylinders. The individual 
cylinders are assumed to be homogeneous and isotropic, just as for option RCTNGL (see ^19.1 8b . 
Let us refer to a single cylinder as the Target Unit Cell (TUC). The TUC is then repeated in the target y- 
and/or z-directions, with periodicities PYDxd and PZDxd, where d is the lattice spacing. To repeat in 
only one direction, set either P YD or PZD to zero. 

This option requires 5 shape parameters: The pertinent line in ddscat . par should read 

SHPARi SHPAR 2 SHPAR 3 SHPAR 4 SHPAR 5 

where SHPARi, SHPAR 2 , SHPAR 3 , SHPAR 4 , SHPAR 5 are numbers: 
SHPARi = cylinder length along axis (in units of d) in units of d 
SHPAR 2 = cylinder diameter/rf 
SHPAR 3 = 1 for cylinder axis || xtf 

= 2 for cylinder axis || yxF 

= 3 for cylinder axis || ztf (see below) 
SHPAR4 = PYD = periodicityAi in ytf direction ( = to suppress repetition) 
SHPAR5 = PZD = periodicityAi in ztf direction ( = to suppress repetition) 

The overall size of the TUC (in terms of numbers of dipoles) is determined by parameters SHPARI 
and SHPAR2. The orientation of a single cylinder is determined by SHPAR3. The periodicity in the TF 
y and z directions is determined by parameters SHPAR4 and SHPAR5. 

The physical size of the TUC is specified by the value of a e s (in physical units, e.g. cm), specified in 
the file ddscat . par, with the usual correspondence d = (47r/3iV) 1 / 3 a c ff, where N is the number of 
dipoles in the TUC. 

With target option CYLNDRPBC, the target becomes a periodic structure, of infinite extent. 

• If NPY > and NPZ > 0, then the target cylindrical TUC repeats in the yxF direction, with 
periodicity NPY x d. 

• If NPY = and NPZ > then the target cylindrical TUC repeats in the ztf direction, with 
periodicity NPZ x d. 

• If NPY > and NPZ > then the target cylindrical TUC repeats in the yxF direction, with 
periodicity NPY x d, and in the ztf direction, with periodicity NPZ x d. 

Target Orientation: The target axes (in the TF) are set to ai — (1, 0, 0)tf a nd &2 — (0, 1, 0)tf- Note 
that ai does not necessarily coincide with the cylinder axis: individual cylinders may have any of 3 
different orientations in the TF. 



Example 1: One could construct a single infinite cylinder with the following two lines in ddscat .par: 
100 1 100 

1.0 100.49 2 1.0 0. 

The first line ensures that there will be enough memory allocated to generate the target. The TUC 
would be a thin circular "slice" containing just one layer of dipoles. The diameter of the circular slice 
would be about 100.49d in extent, so the TUC would have approximately (tt/4) x (100. 49) 2 = 7931 
dipoles (7932 in the actual realization) within a 100 x 1 x 100 "extended target volume". The TUC would 
be oriented with the cylinder axis in the yxF direction (SHPAR3=2) and the structure would repeat in 
the yxF direction with a period of 1.0 x d. SHPAR5=0 means that there will be no repetition in the z 
direction. As noted above, the "target axis" vector ai = xtf- 

Note that SHPARi, SHPAR 2 , SHPAR 4 , and SHPAR 5 need not be integers. However, SHPAR 3 , deter- 
mining the orientation of the cylinders in the TF, can only take on the values 1,2,3. 
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The orientation of the incident radiation relative to the target is, as for all other targets, set by the 
usual orientation angles (3, 6, and $ (see iJT7]above); for example, 9 = would be for radiation incident 
normal to the periodic structure. 

Sample calculation: The subdirectory examples_exp/CYLNDRPBC contains ddscat .par for 
calculating scattering by an infinite cylinder with m = 1.33 + O.Oli for 2irR/\ = 10, where R is the 
cylinder radius. 

20.5 DSKBLYPBC = Target consisting of a periodic array of disks on top of a 
two-layer rectangular slabs. 

This option causes DDSCAT to create a target consisting of a periodic or biperiodic array of Target Unit 
Cells (TUCs), each TUC consisting of a disk of composition 1 resting on top of a rectangular block con- 
sisting of two layers: composition 2 on top and composition 3 below. Materials 1, 2, and 3 are assumed 
to be homogeneous and isotropic. 
This option requires 8 shape parameters: 
The pertinent line in ddscat . par should read 

SHPARi SHPAR.2 SHPAR 3 SHPAR 4 SHPAR 5 SHPAR 6 SHPAR 7 SHPAR 8 
where 

SHPARi = disk thickness in x direction (in Target Frame) in units of d 

SHPAR 2 = (disk diameter)/d 

SHPAR3 = (upper slab thickness)/d 

SHPAR 4 = (lower slab thickness)/c? 

SHPAR5 = (slab extent in yxF direction)/^ 

SHPARg = (slab extent in ztf direction)/<i 

SHPAR7 = period in yxF direction/rf 

SHPARg = period in ztf direction/rf 

The physical size of the TUC is specified by the value of a e g (in physical units, e.g. cm), specified 
in the file ddscat . par. 

The "computational volume" is determined by 

(SHPARi + SHPAR 3 + SHPAR 4 )XSHPAR 5 XSHPAR 6 . 

The lower surface of the cylindrical disk is in the x = plane. The upper surface of the slab is also 
in the x = plane. It is required that SHPAR 2 < min(SHPAR 4 , SHPAR5). 

The Target Frame origin (0,0,0) is located where the symmetry axis of the disk intersects the x = 
plane (the upper surface of the slab, and the lower surface of the disk). 
In the Target Frame, dipoles representing the disk are located at 

x/d = 0.5, 1.5, [int(SHPARi + 0.5) - 0.5] 
and at (y, z) values 

y/d = ±0.5,±1.5,...and 

z/d = ±0.5, ±1.5, ... satisfying 

{y 2 + z 2 ) < (SHPAR 2 /2) 2 d 2 . 
Dipoles representing the rectangular slab are located at (x/d, y/d, z/d) = (j x + 0.5, j y + Ay, j z + 
Az), where j x , j y , and j z are integers. Ay = or 0.5 depending on whether SHPAR5 is even or odd. 
Az = or 0.5 depending on whether SHPARg is even or odd. 

The TUC is repeated in the target y- and z-directions, with periodicities SWKRjxd and SHPARg xd. 

As always, the physical size of the target is fixed by specifying the value of the effective radius 
a ff = (3Vtuc/47t) 1 ^ 3 , where Vttjc is the total volume of solid material in one TUC. For this geometry, 
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the number of dipoles in the target will be approximately 

N = (tt/4) x SHPARi x (SHPAR 2 ) 2 + [SHPAR 3 + SHPAR 4 ] x SHPAR 5 x SHPAR 6 

although the exact number may differ because of the dipoles are required to be located on a rectangular 
lattice. The dipole spacing d in physical units is determined from the specified value of a e ff and the 
number N of dipoles in the target: d = (47r/3iV) 1 ' 3 a ff . 

The target axes (in the TF) are set to a.\ = xtf = (1, 0, 0)tf - i-e., normal to the "slab" - and a.2 = 
yxF = (0, 1, 0)tf- The orientation of the incident radiation relative to the target is, as for all other 
targets, set by the usual orientation angles f3, 0, and $ (see fJTT] above); for example, = would be 
for radiation incident normal to the slab. 

20.6 DSKRCTPBC = Target consisting of homogeneous rectangular brick plus 
a disk, extended in target y and z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of a biperiodic array of Target Unit Cells. 
Each Target Unit Cell (TUC) consists of a disk of composition 1 resting on top of a rectangular block of 
composition 2. Materials 1 and 2 are assumed to be homogeneous and isotropic. 

The cylindrical disk has thickness SHPARi x<i in the x-direction, and diameter SHPAR2XG?. The rect- 
angular block is assumed to have thickness SHPAR.3xd in the x-direction, extent SHPAR 4 x<i in the 
y-direction, and extent SHPAR5XG? in the z-direction. The lower surface of the cylindrical disk is 
in the x = plane. The upper surface of the slab is also in the x = plane. It is required that 
SHPAR 2 <min(SHPAR4, SHPAR 5 ). 

The Target Frame origin (0,0,0) is located where the symmetry axis of the disk intersects the x = 
plane (the upper surface of the slab, and the lower surface of the disk). In the Target Frame, dipoles 
representing the rectangular block are located at (x/d, y/d, z/d) = (j x +0.5, j y + Ay, j z + Az), where 
jx, jy, and j z are integers. Ay = or 0.5 depending on whether SHPAR4 is even or odd. Az = or 0.5 
depending on whether SHPAR5 is even or odd. 
Jx = -[int(SHPAR 3 + 0.5)], -1. 

j y = -[int(0.5 x SHPAR4 - 0.5) + 1] , ... , j z = -[int(SHPAR 4 + 0.5) - 0.5] 

y/d = -[int(0.5 x SHPAR 4 + 0.5) - 0.5], [int(0.5 x SHPAR 4 + 0.5) - 0.5] 

z/d = -[int(0.5 x SHPAR 5 + 0.5) - 0.5], [int(0.5 x SHPAR 5 + 0.5) - 0.5] 
where int(x) is the greatest integer less than or equal to x. Dipoles representing the disk are located at 

x/d = 0.5, 1.5, [int(SHPAR 4 + 0.5) - 0.5] 
and at (y, z) values 

y/d = ±0.5, ±1.5, ... and 

z/d = ±0.5, ±1.5, ... satisfying 

{y 2 + z 2 ) < (SHPAR 5 /2) 2 cP. 
The TUC is repeated in the target y- and z-directions, with periodicities SHPARg x d and SHPAR7 x d. 
As always, the physical size of the target is fixed by specifying the value of the effective radius a c g = 
(3Vtuc/47t) 1 ^ 3 , where Vtuc is the total volume of solid material in one TUC. For this geometry, 
the number of dipoles in the target will be approximately TV = [SHPARi x SHPAR 2 x SHPAR3 + 
(7r/4)((SHPAR 4 ) 2 x SHPAR5)], although the exact number may differ because of the dipoles are re- 
quired to be located on a rectangular lattice. The dipole spacing d in physical units is determined from 
the specified value of a c s and the number N of dipoles in the target: d = (47r/3A r ) 1 / 3 a e ff. This option 
requires 7 shape parameters: 
The pertinent line in ddscat . par should read 

SHPARi SHPAR 2 SHPAR3 SHPAR 4 SHPAR5 SHPAR 6 SHPAR7 
where 

SHPARi = [disk thickness (in xtf direction)]/^ [material 1] 
SHPAR 2 = (disk diameter)/d 

SHPAR3 = (brick thickness in xtf direction)/d [material 2] 
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SHPAR4 = (brick length in yxF direction)/^ 
SHPAR5 = (brick length in ztf direction)/rf 
SHPARg = periodicity in yxF direction/d 
SHPAR7 = periodicity in ztf directionAi 

The overall extent of the TUC (the "computational volume") is determined by parameters ( SHPAR1 + SHPAR4 ) , 
max ( SHPAR 2 , SHPAR 4 , and max ( SHPAR 3 , SHPAR 5 ) . The periodicity in the TF y and z directions is 
determined by parameters SHPAR6 and SHPAR7. 

The physical size of the TUC is specified by the value of a c g (in physical units, e.g. cm), specified in the 
file ddscat . par. 

The target is a periodic structure, of infinite extent in the target y- and z- directions. The target axes 
(in the TF) are set to ai = xtf = (1,0, 0)tf - i-e-, normal to the "slab" - and a.2 = ytf = (0,1, 0)tf- 
The orientation of the incident radiation relative to the target is, as for all other targets, set by the usual 
orientation angles f3, 0, and $ (see CP71 above): for example, 9 = would be for radiation incident 
normal to the slab. 

Sample calculation: Subdirectory examples_exp/DSKRCTPBC contains ddscat . par for cal- 
culating scattering by a 0.0500/im thick SiaN4 slab with a periodic array of Au disks, with period- 
icity 0.0800/im, disk diameter 0.0500/im, and disk height 0.0400/im, for light with wavelength A = 
0.5320/mi. The radiation is incident at an angle of 60° relative to the surface normal. 

20.7 HEXGONPBC = Target consisting of homogeneous hexagonal prism re- 
peated in target y and/or z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of a periodic or biperiodic array of hexagonal 
prisms. The individual prisms are assumed to be homogeneous and isotropic, just as for option RCTNGL 
(see gjyjl. 

Let us refer to a single hexagonal prism as the Target Unit Cell (TUC). The TUC is then repeated in 
the target y- and z-directions, with periodicities PYDxrf and PZDxd, where d is the lattice spacing. To 
repeat in only one direction, set either P YD or P ZD to zero. 

This option requires 5 shape parameters: The pertinent line in ddscat . par should read 

SHPARi SHPAR 2 SHPAR3 SHPAR4 SHPAR5 

whereSHPARi, SHPAR 2 , SHPAR3, SHPAR4, SHPAR5 are numbers: 
SHPARi = prism length along prism axis (in units of d) in units of d 
SHPAR2 = 2xlength of one hexagonal side/d 

SHPAR3 = 1,2,3,4,5 or 6 to specify prism orientation in the TF (see below) 
SHPAR 4 = PYD = periodicity in TF y direction/d 
SHPAR5 = PZD = periodicity in TF z direction/d 

The overall size of the TUC (in terms of numbers of dipoles) is determined by parameters SHPARi, 
SHPAR2, and SHPAR3. The periodicity in the TF y and z directions is determined by parameters SHPAR4 
and SHPAR5. 

The physical size of the TUC is specified by the value of a e ff (in physical units, e.g. cm), specified in 
the file ddscat . par, with the usual correspondence d = (47r/3A^) 1 / 3 a c ff, where N is the number of 
dipoles in the TUC. 

With target option HEXGONPBC, the target becomes a periodic structure, of infinite extent in the 
target y- and z- directions (assuming both NPY and NPZ are nonzero). 

The target axes (in the TF) are set to £4 = xtf = (1, 0, 0)tf - i- e -> normal to the "slab" - and 
a 2 = Ytf = (0, 1,0) T F- 

The individual hexagons may have any of 6 different orientations relative to the slab: Let unit vectors 
h be || to the axis of the hexagonal prism, and let unit vector f be normal to one of the rectangular faces 
of the hexagonal prism. Then 

SHPAR3 = 1 for h || xtf, f || YTF 
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SHPAR3=2 for h 
SHPAR3 = 3 for h 
SHPAR3 = 4 for h 
SHPAR3 = 5 for h 
SHPAR3 = 6 for h 



x TF , f || ztf 
y T F, f || X T F 
y TF , f || z TF 
ztf, f || x T F 
ztf, f || yTF 



For example, one could construct a single infinite hexagonal column with the following line in 

ddscat . par: 



2.0 100.0 3 2.0 0. 



The TUC would be a thin hexagonal "slice" containing two layers of dipoles. The edges of the hexagon 
would be about 50d in extent, so the TUC would have approximately (3V3/2) x 50 2 x 2 = 12990 
dipoles (13024 in the actual realization) within a 90 x 2 x 100 "extended target volume". The TUC 
would be oriented with the hexagonal axis in the yTF direction, with ztf normal to a rectangular faces 
of the prism (SHPAR3=3), and the structure would repeat in the yTF direction with a period of 2 x d 
(SHPAR4=2 . 0). SHPAR5 = means that there will be no repetition in the z T f direction. 

Note that SHPARi, SHPAR 2 , SHPAR 4 , and SHPAR 5 need not be integers. However, SHPAR 5 , deter- 
mining the orientation of the prisms in the TF, can only take on the values 1,2,3,4,5,6. 

Important Note: For technical reasons, P YD and PZD must not be smaller than the "extended" 
target extent in the yTF and ztf directions. When the GPFAFT option is used for the 3 -dimensional 
FFT calculations, the extended target volume always has dimensions/<i = 2 a 3 b 5 c , where a, b, and c are 
nonnegative integers, with (dimension/d) > 1). 

The orientation of the incident radiation relative to the target is, as for all other targets, set by the 
usual orientation angles (3, 9, and $ (see SjT7]above); for example, 9 = would be for radiation incident 
normal to the periodic structure. 



20.8 LYRSLBPBC = Target consisting of layered slab, extended in target y and 
z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of an array of multilayer bricks, layered in 
the xtf direction. The size of each brick in the yTF and ztf direction is specified. Up to 4 layers are 
allowed. 

The bricks are repeated in the yTF and ztf direction with a specified periodicity. If L y = P y and 
L z = P z , then the target consists of a continuous multilayer slab. For this case, it is most economical to 

set Ly/d = L z /d = P y /d = P z /d = 1. 

If P y = 0, then repetition in the yTF direction is suppressed - the target repeats only in the ztf direc- 
tion. 

If P z = 0, then repetition in the ztf direction is suppressed - the target repeats only in the y^p direction. 
The upper surface of the slab is asssume to be located at xtf = 0. The lower surface of the slab is at 

.XTF = ~Lx = - SHPARI X d. 

The multilayer slab geometry is specified with 9 parameters. The pertinent line in ddscat . par should 
read 



SHPARi SHPAR 2 SHPAR 3 SHPAR4 SHPAR5 SHPAR 6 SHPAR7 SHPAR 8 SHPARg 



where 

SHPARi = L x /d = (brick thickness in xtf direction)/^ 
SHPAR2 = L y /d = (brick extent in yTF direction)/d 
SHPAR3 = L z /d= (brick extent in ztf direction)/rf 
SHPAR4 = fraction of the slab with composition 1 
SHPAR5 = fraction of the slab with composition 2 
SHPARg = fraction of the slab with composition 3 
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SHPAR7 = fraction of the slab with composition 4 

SHPARg = P y /d = (periodicity in yxF direction)/<i (0 to suppress repetition in yxF direction) 
SHPARg = P z /d= (periodicity in ztf direction)/^ (0 to suppress repetition in ztf direction) 

For a slab with only one layer, set SHPAR 5 = 0, SHPAR 6 = 0, SHPAR 7 = 0. 
For a slab with only two layers, set SHPAR 6 = and SHPAR7 = 0. 
For a slab with only three layers, set SHPAR7 = 0. 

The user must set NCOMP equal to the number of nonzero thickness layers. 

The number N of dipoles in one TUC is N = nint(SHPARi) x nint(SHPAR 2 ) x nint(SHPAR 3 ). 

The fractions SHPAR 2 , SHPAR3, SHPAR4, and SHPAR 5 must sum to 1. The number of dipoles in 
each of the layers will be integers that are close to nint(SHPARi * SHPAR4), nint(SHPARi * SHPAR 5 ), 
nint(SHPARi * SHPAR 6 ), nint(SHPARi * SHPAR7), 

The physical size of the TUC is specified by the value of a e ff (in physical units, e.g. cm), specified 
in the file ddscat . par. Because of the way a e s is defined (47ra^ ff /3 = Nd 3 ), it should be set to 

a cS = (S/A^N^d = (3/4^ 3 (L x L y L z ) 1 / 3 

With target option LYRSLBPBC, the target becomes a periodic structure, of infinite extent in the 
target y- and z- directions. The target axes (in the TF) are set to ai = (1, 0, 0)tf - i.e., normal to the 
"slab" - and a 2 = (0,1, 0)tf- The orientation of the incident radiation relative to the target is, as for all 
other targets, set by the usual orientation angles j3, 0, and $ (see CP71 above); for example, = would 
be for radiation incident normal to the slab. 

For this option, there are only two allowed scattering directions, corresponding to transmission and 
specular reflection. DDSCAT will calculate both the transmission and reflection coefficients. 
The last two lines in ddscat . par should appear as in the following example ddscat . par file. This 
example is for a slab with two layers: the slab is 26 dipole layers thick; the first layer comprises 76.92% 
of the thickness, the second layer 23.08% of the thickness. The wavelength is 0.532/im, the thickness is 

L x = (4 7 r/3) 1 / 3 7V 2 2/3 a cff = (4tt/3) 1 / 3 (26) 2 / 3 0.009189 = 0.1300/mi. 
The upper layer thickness is 0.2308L X = 0.0300/im 

ddscat . par is set up to calculate a single orientation: in the Lab Frame, the target is rotated through 
an angle = 120°, with $ = 0. In this orientation, the incident radiation is propagating in the 
(—0.5, 0.866, 0) direction in the Target Frame, so that it is impinging on target layer 2 (Au). 

' =========== Parameter file for v7 . 1 . ===================' 

'**** PRELIMINARIES ****' 

' NOTORQ' = CMTORQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 

' PBCGS2 ' = CMDSOL*6 (PBCGS2 , PBCGST, PETRKP) — select solution method 

' GPFAFT ' = CMETHD* 6 (GPFAFT, FFTMKL) 

' GKDLDR' = CALPHA*6 (GKDLDR, LATTDR) 

' NOTBIN' = CBINFLAG (ALLBIN, ORIBIN, NOTBIN) 

Initial Memory Allocation ****' 
26 1 1 = upper bounds on size of TUC 

Target Geometry and Composition 
' LYRSLBPBC ' = CSHAPE*9 shape directive 

26 1 1 0.7692 0.2308 11= shape parameters SHPAR1 - SHPAR9 

2 = NCOMP = number of dielectric materials 

' /u/draine/work/DDA/diel/Eagle_2 00 ' = refractive index 1 

' /u/draine/work/DDA/diel/Au_evap' = refractive index 2 

'**** Error Tolerance ****' 

1.00e-5 = TOL = MAX ALLOWED (NORM OF | G>=AC I E>-ACA | X> ) / (NORM OF AC|E>) 

'**** interaction cutoff parameter for PBC calculations ****' 

5.00e-3 = GAMMA (le-2 is normal, 3e-3 for greater accuracy) 

'**** Angular resolution for calculation of <cos>, etc. ****' 

0.5 = ETASCA (number of angles is proportional to [ (3+x) /ETASCA] "2 ) 

Wavelengths (micron) ****' 
0.5320 0.5320 1 ' INV = wavelengths ( first , last , how many , how=LIN, INV, LOG) 
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'**** Effective Radii (micron) **** ' 

0.009189 0.009189 1 ' LIN' = eff. radii (first, last, how many , how=LIN, INV, LOG) 
'**** Define Incident Polarizations ****' 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 
2 = IORTH (=1 to do only pol . state eOl; =2 to also do orth. pol . state) 

Specify which output files to write 
1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient. 

I = IWRPOL (=0 to suppress, =1 to write ".pol" file for each (BETA, THETA) 
'**** Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 

120. 120. 1 = THETMI, THETMX, NTHETA (theta=angle between al and k) 

0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 

'**** Specify first I WAV, IRAD, IORI (normally 0) ****' 

= first I WAV, first IRAD, first IORI (0 to begin fresh) 
'**** Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 

II 12 21 22 31 41 = indices ij of elements to print 
'**** Specify Scattered Directions ****' 

' TFRAME ' = CMDFRM (LFRAME, T FRAME for Lab Frame or Target Frame) 

1 = number of scattering orders 
0. 0. = (M,N) for scattering 



20.9 RCTGL_PBC = Target consisting of homogeneous rectangular brick, ex- 
tended in target y and z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of a biperiodic array of rectangular bricks. The 
bricks are assumed to be homogeneous and isotropic, just as for option RCTNGL (see ^19. 18b - 
Let us refer to a single rectangular brick as the Target Unit Cell (TUC). The TUC is then repeated 
in the ?/tf- and ZTF-directions, with periodicities PYAEFFxa e ff and PZAEFFxa e ff, where a c g = 
(3VTUC/47T) 1 / 3 , where Vtuc is the total volume of solid material in one TUC. This option requires 
5 shape parameters: 

The pertinent line in ddscat . par should read 
SHPARi SHPAR.2 SHPAR 3 SHPAR 4 SHPAR 5 
where 

SHPARi = (brick thickness)/d in the xtf direction 
SHPAR2 = (brick thickness)/d in the i/tf direction 
SHPAR3 = (brick thickness)/^ in the ztf direction 
SHPAR4 = periodicity/d in the i/tf direction 
SHPAR5 = periodicity/d in the ztf direction 

The overall size of the TUC (in terms of numbers of dipoles) is determined by parameters SHPARi, 
SHPAR2, and SHPAR3. The periodicity in the yxF and ztf directions is determined by parameters 

SHPAR 4 and SHPAR 5 . 

The physical size of the TUC is specified by the value of a c g (in physical units, e.g. cm), specified in the 
file ddscat . par. 

With target option RCTGL_PBC, the target becomes a periodic structure, of infinite extent in the 
target y- and z- directions. The target axes (in the TF) are set to ai = (1, 0, 0)tf - i.e., normal to the 
"slab" - and &2 = (0,1, 0)tf- The orientation of the incident radiation relative to the target is, as for all 
other targets, set by the usual orientation angles j3, 0, and $ (see iJT7]above); for example, = would 
be for radiation incident normal to the slab. 

Sample Calculation: Subdirectory examples_exp/RCTGL_PBC contains ddscat .par for scat- 
tering by an infinite slab, constituted from 16 x 1 x 1 dipole TUCs. 
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20.10 RECRECPBC = Rectangular solid resting on top of another rectangular 
solid, repeated periodically in target y and z directions using periodic boundary 
conditions 

The TUC consists of a single rectangular "brick", of material 1, resting on top of a second rectangular 
brick, of material 2. The centroids of the two bricks along a line in the xtf direction. The bricks are as- 
sumed to be homogeneous, and materials 1 and 2 are assumed to be isotropic. The TUC is then repeated 
in the utf- and z tf -directions, with periodicities SHPAR4 x d and SHPAR 5 x d. is the total volume of 
solid material in one TUC. This option requires 5 shape parameters: 
The pertinent line in ddscat . par should read 



SHPARi SHPAR 2 SHPAR3 SHPAR4 SHPAR5 



where 

SHPARi = (upper brick thickness)/<i in the xtf direction 
SHPAR2 = (upper brick thickness)/rf in the yxF direction 
SHPAR 3 = (upper brick thickness)/rf in the ztf direction 
SHPAR4 = (lower brick thickness)/rf in the xtf direction 
SHPAR5 = (lower brick thickness)/d in the yxF direction 
SHPARg = (lower brick thickness)/rf in the ztf direction 
SHPAR7 = periodicity/cJ in the yxF direction 
SHPARs = periodicity Id in the ztf direction 



The actual numbers of dipoles Nix, N\ y , N\ z , along each dimension of the upper brick, and N 2x , 
N 2y , N 2z along each dimension of the lower brick, just be integers. Usually, N\ x = nint(SHPARi), 
Ni y = nint(SHPAR 2 ), N Xz = nint(SHPAR 3 ), N 2x = nint(SHPAR 4 ), N 2y = nint(SHPAR 5 ), N 2z = 
nint(SHPAR6), where ninths) is the integer nearest to x, but under some circumstances Ni x , N\ y , Ni z , 
N 2x , N 2y , N 2z might be larger or smaller by 1 unit. 

The overall size of the TUC (in terms of numbers of dipoles) is determined by parameters SHPARi 
- SHPAR 6 : 

N = (N lx x Ni y xN lz ) + {N 2x x N 2y x N 2z ) (34) 

The periodicity in the yxF and ztf directions is determined by parameters SHPAR7 and SHPARg. 
The physical size of the TUC is specified by the value of a e g (in physical units, e.g. cm), specified in the 
file ddscat . par: 

d = (47r/3A0 1/3 a off (35) 

The target is a periodic structure, of infinite extent in the yxF and ztf directions. The target axes are 
set to a.i = xtf - i.e., normal to the "slab" - and k 2 = yTF- The orientation of the incident radiation 
relative to the target is, as for all other targets, set by the usual orientation angles /3, 8, and $ (see SjTT] 
above) specifying the orientation of the target axes ai and k 2 relative to the direction of incidence; for 
example, 6 = would be for radiation incident normal to the slab. 

The scattering directions are specified by specifying the diffraction order (M, N); for each diffrac- 
tion order one transmitted wave direction and one reflected wave direction will be calculated, with the 
dimensionless 4x4 scattering matrix S^ 2 " 1 ) calculated for each scattering direction. At large distances 
from the infinite slab, the scattered Stokes vector in the (M,N) diffraction order is 

4 

I sca ,i (M, N) = Y, S^hnj (36) 

where Ii n j is the incident Stokes vector. See Draine & Flatau (2008) for interpretation of the Sij as 
transmission and reflection efficiencies. 
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20. 1 1 SLBHOLPB C = Target consisting of a periodic array of rectangular blocks, 
each containing a cylindrical hole 

Individual blocks have extent (a,b,c) in the (xtf, Ytf, ztf) directions, and the cylindrical hole has 
radius r. The period in the yxF-direction is P y , and the period in the zxF-direction is P z . 
The pertinent line in ddscat . par should consist of 

SHPARi SHPAR 2 SHPAR 3 SHPAR 4 SHPAR 5 SHPAR 6 

where SHPARi = a/d(d is the interdipole spacing) 

SHPAR 2 = b/a 

SHPAR 3 = c/a 

SHPAR 4 = r/a 

SHPAR 5 = P y /d 

SHPARg = P z /d. 

Ideally, a/d =SHPARi, b/d =SHPAR 2 x SHPARi, and c/d =SHPAR 3 x SHPARi will be integers (so 
that the cubic lattice can accurately approximate the desired target). If P y = and P z > the target is 
periodic in the ztf -direction only. 

If P y > and P z = the target is periodic in the yxF -direction only. 

If P y > it is required that P y > b, and if P z > it is required that P z > c, so that the blocks do not 
overlap. 

With P y = b and P z = c, the blocks are juxtaposed to form a periodic array of cylindrical holes in a 
solid slab. 

Example: ddscat_SLBHOLPBC .par is a sample ddscat .par for a periodic array of cylindrical 
holes in a slab of thickness a, with holes of radius r, and period P y and Pd 

20.12 SPHRN_PBC = Target consisting of group of N spheres, extended in target 
y and z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of a periodic array of N— sphere structures, 
where one N— sphere structure consists of N spheres, just as for target option NANSPH (see i)19.23| >. 
Each sphere can be of arbitary composition, and can be anisotropic if desired. Information for the de- 
scription of one TV-sphere structure is supplied via an external file, just as for target option NANSPH - 
see SfljOS). 

Let us refer to a single TV-sphere structure as the Target Unit Cell (TUC). The TUC is then repeated 
in the i/tf- and ZTF-directions, with periodicities PYAEFF xa c ff and PZAEFFxa c ff, where a c s = 
(3VTUC/47T) 1 / 3 , where Vtuc is the total volume of solid material in one TUC. This option requires 
3 shape parameters: 

D I AMX = maximum extent of target in the target frame x direction/d 
PYAEFF = periodicity in target y direction/a c ff 
PZAEFF = periodicity in target z direction/a c ff. 
The pertinent line in ddscat. par should read 

SHPARi SHPAR 2 SHPAR3 'filename' (quotes must be used) 

where 

SHPARi = target diameter in x direction (in Target Frame) in units of d 

SHPAR 2 = PYAEFF 
SHPAR 3 = PZAEFF. 

filename is the name of the file specifying the locations and relative sizes of the spheres. 

The overall size of the TUC (in terms of numbers of dipoles) is determined by parameter SHPARi, 
which is the extent of the multisphere target in the x-direction, in units of the lattice spacing d. The 
physical size of the TUC is specified by the value of a B s (in physical units, e.g. cm), specified in the file 

ddscat . par. 
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The location of the spheres in the TUC, and their composition, is specified in file 'filename'. Please 
consult i jl9.23l above for detailed information concerning the information in this file, and its arrangement. 

Note that while the spheres can be anisotropic and of differing composition, they can of course also 
be isotropic and of a single composition, in which case the relevant lines in the file filename ' should 
read 

X! yi zi ri 1 1 1 
x 2 2/2 z 2 r 2 1 1 1 
x 3 y 3 z 3 r 3 1 1 1 

i.e., every sphere has isotropic composition ICOMP = l and the three dielectric function orientation angles 
are set to zero. 

When the user uses target option SPHRN_PBC, the target now becomes a periodic structure, of infinite 
extent in the target y- and z- directions. The target axis ai = xlf = (1, 0, 0)tf - i.e., normal to the 
"slab" - and target axis & 2 = Ylf = (0,1, 0)tf- The orientation of the incident radiation relative to the 
target is, as for all other targets, set by the usual orientation angles f3, 0, and <j> (see fllTl above). The 
scattering directions are, just as for other targets, determined by the scattering angles S , <f> s (see SfJTJ 
below). 

The scattering problem for this infinite structure, assumed to be illuminated by an incident monochro- 
matic plane wave, is essentially solved "exactly", in the sense that the electric polarization of each of the 
constituent dipoles is due to the electric field produced by the incident plane wave plus all of the other 
dipoles in the infinite target. 

However, the assumed target will, of course, act as a perfect diffraction grating if the scattered radia- 
tion is calculated as the coherent sum of all the oscillating dipoles in this periodic structure: the far-field 
scattered intensity would be zero in all directions except those where the Bragg scattering condition is 
satisfied, and in those directions the far-field scattering intensity would be infinite. 

To suppress this singular behavior, we calculate the far-field scattered intensity as though the sepa- 
rate TUCs scatter incoherently. The scattering efficiency Qsca and the absorption efficiency Q a b s are 
defined to be the scattering and absorption cross section per TUC, divided by tto^ s , where a e g = 
(3Vtuc /47T) 1 / 3 , where Vtuc is the volume of solid material per TUC. 

Note: the user is allowed to set the target periodicity in the target j/tf (or ztf) direction to values that 
could be smaller than the total extent of one TUC in the target j/tf (or ztf) direction. This is physically 
allowable, provided that the spheres from one TUC do not overlap with the spheres from neighboring 
TUCs. Note that DDSCAT does not check for such overlap. 

20.13 TRILYRPBC = Three stacked rectangular blocks, repeated periodically 

The target unit cell (TUC) consists of a stack of 3 rectangular blocks with centers on the xtf axis. The 
TUC is repeated in either the yxF direction, the ztf direction, or both. A total of 1 1 shape parameters 
must be specified: 

SHPARi = x-thickness of upper layer/d [material 1] 

SHPAR 2 = y-width/d of upper layer 

SHPAR 3 = z-width/rf of upper layer 

SHPAR4 = x-thicknessAi of middle layer [material 2] 

SHPAR 5 = y-width/d of middle layer 

SHPAR 6 = z-width/rf of middle layer 

SHPAR7 = x-width/<i of lower layer 

SHPAR 8 = y-width/d of lower layer 

SHPARg = z-width/rf of lower layer 

SHPAR10 = period/tf in y direction 

SHPARn= period/d in z direction 
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21 Scattering Directions 
21.1 Isolated Finite Targets 

DDSCAT calculates scattering in selected directions, and elements of the scattering matrix are reported 
in the output files wxxxryyykzzz . sea . The scattering direction is specified through angles 8 S and <f) s 
(not to be confused with the angles 9 and $ which specify the orientation of the target relative to the 
incident radiation!). 

For isolated finite targets (i.e., PBC not employed) there are two options for specifying the scatter- 
ing direction, with the option determined by the value of the string CMDFRM read from the input file 

ddscat . par. 

1. If the user specifies CMDFRM=' LFRAME', then the angles 9, <f> input from ddscat .par are 
understood to specify the scattering directions relative to the Lab Frame (the frame where the 
incident beam is in the x— direction). 

When CMDFRM=' LFRAME' , the angle 9 is simply the scattering angle 9 S : the angle between the 
incident beam (in direction xlf) and the scattered beam (9 S = for forward scattering, S = 180° 
for backscattering). 

The angle <fi specifies the orientation of the "scattering plane" relative to the xlf — Ylf plane. 
When = the scattering plane is assumed to coincide with the xlf — Ylf plane. When <f> = 90° 
the scattering plane is assumed to coincide with the xlf — zlf plane. Within the scattering plane 
the scattering directions are specified by < 9 < 180°. Thus: 



2. If the user specifies CMDFRM=' TFRAME ' , then the angles 0, (f> input from ddscat .par are 
understood to specify the scattering directions n s relative to the Target Frame (the frame defined 
by target axes xtf, Ytf, ztf)- is the angle between n s and xtf, and <\> is the angle between the 
n s — xtf plane and the xtf — Ytf plane. Thus: 



Scattering directions for which the scattering properties are to be calculated are set in the parameter 
file ddscat . par by specifying one or more scattering planes (determined by the value of 4> s ) and for 
each scattering plane, the number and range of 9 S values. The only limitation is that the number of 
scattering directions not exceed the parameter MXSCA in DDSCAT . f (in the code as distributed it is set 

to MXSCA=1000). 

21.2 Scattering Directions for Targets that are Periodic in 1 Dimension 

For targets that are periodic, scattering is only allowed in certain directions (see Draine & Flatau 2008). 
If the user has chosen a PBC target (e.g, CYLNDRPBC, HEXGONPBC, or RCTGL_PBC), SPHRN.PBC 
but has set one of the periodicities to zero, then the target is periodic in only one dimension - e.g., 
CYLNDRPBC could be used to construct a single infinite cylinder. 

In this case, the scattering directions are specified by giving an integral diffraction order M = 
0,±1,±2, ... and one angle, the azimuthal angle ( around the target repetition axis. The diffraction 
order M determines the projection of k s onto the repetition direction. For a given order M, the scatter- 
ing angles with £ = —> 2n form a cone around the repetition direction. 

For example, if PYD > (target repeating in the y^p direction), then M determines the value of 
k sy = k s • ytf, where yxF is the unit vector in the Target Frame y— direction: 



n s = xlf cos 9 + ylf sin 9 cos <f> + zlf sin 9 sin <fi 



(37) 



n s = xtf cos 9 + yxF sin 9 cos <p + ztf sin 9 sin 



(38) 



k sy = ko y + 2t:M/L. 



v 



(39) 



where k 0y = k • Ytfi where k is the incident k vector. Note that the diffraction order M must satisfy 
the condition 

(k 0y - k )(L y /27r) < M < (fc - k y)(L y /2Tr) (40) 



SCATTERING DIRECTIONS 



where L y = PYD x d is the periodicity along the y axis in the Target Frame. M = is always an allowed 
diffraction order. 

The azimuthal angle ( defines a right-handed rotation of the scattering direction around the target 
repetition axis. Thus for a target repetition axis ytf, 

k sx = k ± cos( , (41) 
k sz = fcxsinC , (42) 

where k± = (k 2 — k 2 ^) 1 ^ 2 , with k sy = ko y + 2nM/L y . For a target with repetition axis ztf, 

k sx = k±cos( , (43) 
k sy = fcxsinC , (44) 

where k ± = (k 2 - k 2 ,) 1 / 2 , k sz = k 0z + 2ttM/L z . 

The user selects a diffraction order M and the azimuthal angles £ to be used for that M via one line 
in ddscat . par. An example would be to use CYLNDRPBC to construct an infinite cylinder with the 
cylinder direction in the yTF direction: 

' ============== Parameter file for v7 . 1 . =================== ' 

* **** PRELIMINARIES ****' 

' NOTORQ' = CMT0RQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 

' PBCGS2 ' = CMDS0L*6 (PBCGS2, PBCGST, PETRKP) — select solution method 

' GPFAFT' = CMETHD*6 (GPFAFT, FFTMKL) 

'GKDLDR' = CALPHA*6 (GKDLDR, LATTDR) 

'NOTBIN' = CBINFLAG (NOTBIN, ORIBIN, ALLBIN) 

' **** Initial Memory Allocation ***** 

50 50 50 = size allowance for target generation 

r **** Target Geometry and Composition ****' 

'CYLNDRPBC' = CSHAPE*9 shape directive 

2.0 20. 2 2.0 0. = shape parameters SHPAR1, SHPAR2, SHPAR3, SHPAR4, SHPAR5, 

1 = NCOMP = number of dielectric materials 

' /u/draine/work/DDA/diel/ml . 33_0 . 01' = name of file containing dielectric function 
' **** Error Tolerance ****' 

1.00e-5 = TOL = MAX ALLOWED (NORM OF | G>=AC I E>-ACA | X> ) / (NORM OF AC|E>) 

' **** Interaction cutoff parameter for PBC calculations ***** 

5.00e-3 = GAMMA (le-2 is normal, 3e-3 for greater accuracy) 

' **** Angular resolution for calculation of <cos>, etc. ***** 

0.5 = ETASCA (number of angles is proportional to [ ( 3+x) /ETASCA] " 2 ) 

' **** Wavelengths (micron) ***** 

6.283185 6.283185 1 ' INV = wavelengths (first, last, how many, how=LIN, INV, LOG) 
' **** Effective Radii (micron) **** ' 

2 2 1 'LIN' = eff. radii (first, last, how many, how=LIN, INV, LOG) 
' **** Define Incident Polarizations ***** 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = IORTH (=1 to do only pol . state eOl; =2 to also do orth. pol . state) 

1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient.) 

I = IWRPOL (=0 to suppress, =1 to write ".pol" file for each (BETA, THETA) ) 

* **** Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 

60. 60. 1 = THETMI, THETMX, NTHETA (theta=angle between al and k) 

0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 

'**** Specify first IWAV, IRAD, IORI (normally 0) ****' 

= first IWAV, first IRAD, first IORI (0 to begin fresh) 

* **** Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 

II 12 21 22 31 41 = indices ij of elements to print 
' **** Specify Scattered Directions ***** 

' TFRAME' = CMDFRM (LFRAME, T FRAME for Lab Frame or Target Frame) 

1 = number of scattering cones 

0. 0. 360. 30 = M, zeta_min, zeta_max, dzeta (in degrees) 

In this example, a single diffraction order M = is selected, and £ is to run from f m j n = to 
(max = 360° in increments of 6( = 30°. 

There may be additional lines, one per diffraction order. Remember, however, that every diffraction 
order must satisfy eq. d40b . 
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21.3 Scattering Directions for Targets for Doubly-Periodic Targets 

If the user has specified nonzero periodicity in both the y and z directions, then the scattering direc- 
tions are specified by two integers - the diffraction orders M, N for the y, z directions. The scattering 
directions are 

, k n „ ~ w /, - 27rAA a / „ 2ttN\ „ 

k s = ± — (k ■ x T f) x T f + k ■ ytf H z — ytf + k ■ z T f + — z T f (45) 



(46) 



^0 \ Ly J V L z J 

21 V2 



2ttM\ 2 / „ 2tt^\ 2 



where the + sign gives transmission, and the — sign gives reflection. The integers M and N must 
together satisfy the inequality 

(k 0y + 2-KM/Ly) 2 + {k 0z + 2nN/L z ) 2 < k 2 (47) 

which, for small values of L y and L z , may limit the scattering to only (M,N) = (0,0). [Of course, 
(0, 0) is always allowed]. For each (M, TV) specified in ddscat . par, DDSCAT 7.1 will calculate the 
generalized Mueller matrix s\ 2d ^ (M, N) - see §At large distances from the infinite slab, the scattered 
Stokes vector in the (M, N) diffraction order is 

4 

I sca ^ (M, N)=J2 Sg d) (M, N)I in>j (48) 

3=1 

where Ij n j is the incident Stokes vector, and S 1 ^ (M, N) is the generalization of the 4x4 Miiller 
scattering matrix to targets that are periodic in 2-directions (Draine & Flatau 2008). There are distinct 
Sjj (M, N) for transmission and for reflection, corresponding to the ± in eq. d45l l. 

Draine & Flatau (2008, eq. 69-71) show how the S^ d '(M, N) are easily related to familiar "trans- 
mission coefficients" and "reflection coefficients". 
Here is a sample ddscat . par file as an example: 

' ============== Parameter file for v7 . 1 . ' 

'**** PRELIMINARIES 

' NOTORQ' = CMTORQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 
' PBCGS2 ' = CMDSOL*6 (PBCGS2, PBCGST, PETRKP) — select solution method 
'GPFAFT' = CMETHD*6 (GPFAFT, FFTMKL) 
'GKDLDR' = CALPHA*6 (GKDLDR, LATTDR) 
'NOTBIN' = CBINFLAG (NOTBIN, ORIBIN, ALLBIN) 

Initial Memory Allocation 
50 50 50 = dimensioning allowance for target generation 
'**** Target Geometry and Composition 
'RCTGL_PBC = CSHAPE*9 shape dirctive 

20.0 1.0 1.0 1.0 1.0 = shape parameters PARI, PAR2 , PAR3 , . . . 

1 = NCOMP = number of dielectric materials 

' /u/draine/work/DDA/diel/ml . 33_0 . 1 ' = name of file containing dielectric function 
'**** CONJUGATE GRADIENT DEFINITIONS 

= INIT (TO BEGIN WITH |X0> = 0) 

1.0 0e-5 = TOL = MAX ALLOWED (NORM OF | G>=AC I E>-ACA | X> ) / (NORM OF AC|E>) 

'**** Interaction cutoff parameter for PBC calculations ****' 

5.00e-3 = GAMMA (le-2 is normal, 3e-3 for greater accuracy) 

'**** Angular resolution for calculation of <cos>, etc. ****' 

0.5 = ETASCA (number of angles is proportional to [ (3+x) /ETASCA] ~2 ) 

'**** Wavelengths (micron) ****' 

6.283185 6.283185 1 ' INV = wavelengths (first, last, how many, how=LIN, INV, LOG) 
'**** Effective Radii (micron) ***+ ' 

2 2 1 'LIN' = eff. radii (first, last, how many, how=LIN, INV, LOG) 
' **** Define Incident Polarizations ****' 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = IORTH (=1 to do only pol . state eOl; =2 to also do orth. pol . state) 

1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient.) 
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I = IWRPOL (=0 to suppress, =1 to write ".pol" file for each (BETA, THETA) ) 
'**** Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 
60. 60. 1 = THETMI, THETMX, NTHETA (theta=angle between al and k) 
0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 
Specify first I WAV, IRAD, IORI (normally 0) 

= first I WAV, first IRAD, first IORI (0 to begin fresh) 
' **** Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 

II 12 21 22 31 41 = indices ij of elements to print 
' **** Specify Scattered Directions ****' 

'TFRAME' = CMDFRM (LFRAME, T FRAME for Lab Frame or Target Frame) 

1 = number of scattering orders 
= M, N (diffraction orders) 

22 Incident Polarization State 

Recall that the "Lab Frame" is defined such that the incident radiation is propagating along the xlf 
axis. DDSCAT allows the user to specify a general elliptical polarization for the incident radiation, by 
specifying the (complex) polarization vector eoi. The orthonormal polarization state £02 = xlf x 
is generated automatically if ddscat .par specifies IORTH=2. 

For incident linear polarization, one can simply set eoi = y by specifying 
(0,0) (1,0) (0,0) 

in ddscat . par; then 602 = z. For polarization mode eoi to correspond to right-handed circular po- 
larization, set eoi = (y + iz)/y / 2 by specifying (0,0) (1,0) (0,1) in ddscat . par (DDSCAT 
automatically takes care of the normalization of eoi); then 602 = (if + z)/y/2, corresponding to left- 
handed circular polarization. 

23 Averaging over Scattering Directions: g(l) = (cos 9 S ), etc. 
23.1 Angular Averaging 

An example of scattering by a nonspherical target is shown in Fig. [SJ showing the scattering for a tilted 
cube. Results are shown for four different scattering planes. 

DDSCAT automatically carries out numerical integration of various scattering properties, including 

• (co&0 s ); 

• (cos 2 9 S ); 

• g = (cos 6* s )xlf + (sin S cos 4> s )y + (sin S sin <fi s )z (see §151 : 

• Qr, provided option DOTORQ is specified (see Ci"5l l. 

The angular averages are accomplished by evaluating the scattered intensity for selected scattering di- 
rections (6 S , (f) s ), and taking the appropriately weighted sum. Suppose that we have Ng different values 



For a given j, the values of <fij,k are assumed to be uniformly spaced: <fij,k+i — 4>j,k = 2ir/N < / > (j). The 
angular average of a quantity f(6 s , 4> s ) is approximated by 




(49) 



(50) 




(51) 
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[cos^-j.) - cos(fy +1 )] 
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Ng — 1 



(52) 
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Figure 8: Scattered intensity for incident unpolarized light for a cube with m = 1.5 + O.OOlz and x — 2na c s:/X = 6.5 (i.e., 
d/X = 1.6676, where d is the length of a side). The cube is tilted with respect to the incident radiation, with <d — 30°, and rotated 
by /3 — 15° around its axis to break reflection symmetry. The Mueller matrix element Su is shown for 4 different scattering planes. 
The strong forward scattering lobe is evident. It is also seen that the scattered intensity is a strong function of scattering angle (f> s as 
well as 9 a - at a given value of 9 S (e.g., 9 S = 120°), the scattered intensity can vary by orders of magnitude as <j> s changes by 90°. 
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For a sufficiently large number of scattering directions 

7V 8ca = ^7V^(j) 

3=1 



(53) 
(54) 

(55) 



the sum ( l5TT l approaches the desired integral, but the calculations can be a significant cpu-time burden, 
so efficiency is an important consideration. 



23.2 Selection of Scattering Angles 9 S , d> s 

Beginning with DDSCAT 6.1, an improved approach is taken to evaluation of the angular average. Since 
targets with large values of x = 27ra c ff / A in general have strong forward scattering, it is important to 
obtain good sampling of the forward scattering direction. To implement a preferential sampling of the 
forward scattering directions, the scattering directions (6, <f>) are chosen so that 6 values correspond 
to equal intervals in a monotonically increasing function s(9). The function s(6) is chosen to have a 
negative second derivative d 2 s/d9 2 < so that the density of scattering directions will be higher for 
small values of 9. DDSCAT takes 

s ^ = 6+ nTT ■ (56) 
V + fro 

This provides increased resolution (i.e., increased ds/d9) for 9 < 6q. We want this for the forward 
scattering lobe, so we take 

= t^- . (57) 
1 + x 
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Figure 9: Errors in (cos 6) and (cos 2 9) calculated for a N — 17904 dipole pseudosphere with m = 1.33 + O.Oli, as functions 
of x — 2ira/\. Results are shown for different values of the parameter 77. 
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(58) 
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Thus for small values of x, max[A0] = 30°?7, and for x ^> 1, max[A0] — > 90°77/a;. For a sphere, 
minima in the scattering pattern are separated by ~ 180°/a;, so 77 = 1 would be expected to marginally 
resolve structure in the scattering function. Smaller values of r\ will obviously lead to improved sampling 
of the scattering function, and more accurate angular averages. 

The scattering angles 9j used for the angular averaging are then given by 



where 



(gj -1-0 O )+ [(1 + 9 - Sj ) 2 +W oSj 
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Ng-1 
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(61) 



(62) 



For each 9j, we must choose values of (j>. For 9\ = and 9n = ir only a single value of <f> is needed 
(the scattering is independent of <fi in these two directions). For < 6j < % we use 



N# = max{3,nint [An sm(9j)/ (9 J+1 - ^-1)]} 



(63) 



where nint = nearest integer. This provides sampling in <f> consistent with the sampling in 9. 
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Figure 10: Same as Fig.[9l but for a TV = 32768 dipole cube with m = 1.5 + O.OOli 





1 1 1 


error[(cos0)] 


















_ 3.~=0.125 







23.3 Accuracy of Angular Averaging as a Function of rj 

Figure [9] shows the absolute errors in (cos 8) and (cos 2 6) calculated for a sphere with refractive index 
m = 1.33 + O.Olz using the above prescription for choosing scattering angles. The error is shown as a 
function of scattering parameter x. We see that accuracies of order 0.01 are attained with r\ = 1, and that 
the above prescription provides an accuracy which is approximately independent of x. We recommend 
using values of r\ < 1 unless accuracy in the angular averages is not important. 



24 Scattering by Finite Targets: The Mueller Matrix 
24.1 Two Orthogonal Incident Polarizations (I0RTH=2) 

Throughout the following discussion, x, y, z are unit vector defining the Lab Frame (thus ko = x). 

DDSCAT internally computes the scattering properties of the dipole array in terms of a complex 
scattering matrix f m i(8s, </>s) (Draine 1988), where index 1 = 1,2 denotes the incident polarization state, 
m = 1,2 denotes the scattered polarization state, and 8 s ,<j> s specify the scattering direction. Normally 
DDSCAT is used with IORTH=2 in ddscat .par, so that the scattering problem will be solved for 
both incident polarization states (7 = 1 and 2); in this subsection it will be assumed that this is the case. 

Incident polarization states I = 1, 2 corresp ond t o polarization states e i, e 02 ; recall that polarization 
state e i is user-specified, and e 02 = x x 1( l n Scattered polarization state m = 1 corresponds to 
linear polarization of the scattered wave parallel to the scattering plane (ei = = 6 S ) and m = 2 
corresponds to linear polarization perpendicular to the scattering plane (in the +(j> s direction). The 
scattering matrix f m i was defined (Draine 1988) so that the scattered electric field E s is related to the 
incident electric field Ej(0) at the origin (where the target is assumed to be located) by 



E s - 

E. s • 



exp(ik s • r) 
kr 



hi 



fl2 
$22 



E,(0) 
Ei(0) 



-oi 

i iY2 



(64) 



1 Draine (1988) adopted the convention eo2 = x x eoi. The customary definition of ejx is e;x = &a\ x k. Thus, if eoi 
then eo2 = —etx ■ 

" A frequent choice is eoi = Ylf, with eo2 = xlf x yLF = zlf- 
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The 2x2 complex scattering amplitude matrix (with elements Si, S 2 , S3, and S4) is defined so that (see 
Bohren & Huffman 1983) 



E s ■ §, \ = cxp(ik s ■ r) / S 2 S 3 \ ( Ei(0) • e 4 || 
-E s • L ) -ikr 1 S 4 Si J \ Ei(0) • e l± _ 



(65) 



where e^ii , e;j_ are (real) unit vectors for incident polarization parallel and perpendicular to the scattering 
plane (with the customary definition of e,j_ = e^ii x xlf)- 



From ( 164165b we may write 



S 2 S 3 \ / Ei(0) • e iy \ / hi fi 2 \(Ei(0)-e* 01 

Si Si Jl Ei(0) -ea J V -/21 -/aa 7 I ^(0) • 



(66) 



Let 



a 


e 01 


■yLF , 


(67) 


b = 


^01 


' ZLF , 


(68) 


c = 


^02 


■yLF , 


(69) 


d = 


^02 


' ZLF ■ 


(70) 



Note that since §01, 602 could be complex (i.e., elliptical polarization), the quantities a, b, c, d are com- 
plex. Then 



-01 
-02 



a b \ f y LF 
c d ) V z LF 



and eq. d66l l can be written 



S 2 S 3 \f E,-(0) - e«,| \ f -fu -fi2 \f a b\f E,(0) • y LF 
S 4 S 1 J\ Ei(0)-ecL 7 V /21 /aa A c d / V E *(0)-z LF 

The incident polarization states ejii and e,j_ are related to xlf, yLF, zlf by 

yLF \ _ ( cos (j) s sxa<f> 8 \ ( e 4 || 
zlf / V siring - cos <f> s I \ e i± 



(71) 



(72) 



(73) 



substituting d73l ) into ( 1721 we obtain 

S 2 S 3 \ ( Ei(0) • e,|| \ . / -/ u -/ 12 \( a b\( cos<j> s sin<j> s \ / E*(0) • e,,, 
S 4 5i A E >(°) ' ^ J \ hi /22 J\ c d J { sin^ -cosfa J \ Ei(0) • e u 

(74) 

Eq. d74l i must be true for all Ej(0), so we obtain an expression for the complex scattering amplitude 
matrix in terms of the f m i : 



( s 2 


s 3 " 




\ s 4 


Si j 


M 



-fix -/12 \ f a b \ ( cos0 s sint/), 
/21 f 22 J \ c d ) \ sin^ s -cos< 



(75) 



This provides the 4 equations used in subroutine GETMUELLER to compute the scattering amplitude 
matrix elements: 

51 = -i [f 2 i(b cos 4> s - a sin 4> s ) + f 22 (d cos 4> s - c sin 4> s )] , (76) 

5 2 = -i [f ii(a cos 4> s + b sin <fi s ) + fi 2 (c cos 4> s + d siii 4> s )] , (77) 

53 = i [/n(6cos0 s - asin0 s ) + /i 2 (dcos0 s - csin0 s )] , (78) 

54 = i [f 2 i(a cos (j> s +b sin (j> s ) + f 22 (c cos c/) s + d sin (j> s )} . (79) 
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24.2 Stokes Parameters 

It is both convenient and customary to characterize both incident and scattered radiation by 4 "Stokes 
parameters" - the elements of the "Stokes vector". There are different conventions in the literature; we 
adhere to the definitions of the Stokes vector (I,Q,U,V) adopted in the excellent treatise by Bohren & 
Huffman (1983), to which the reader is referred for further detail. Here are some examples of Stokes 
vectors (I, Q, U, V) = (1,Q/I, U/I, V/I)I: 

• (1,0,0, 0)7: unpolarized light (with intensity I); 

• (1, 1, 0, 0)1 : 100% linearly polarized with E parallel to the scattering plane; 

• (1,-1,0, 0)1 : 100% linearly polarized with E perpendicular to the scattering plane; 

• (1, 0, 1, 0)1 : 100% linearly polarized with E at +45° relative to the scattering plane; 

• (1, 0, —1, 0)1 : 100% linearly polarized with E at -45° relative to the scattering plane; 

• (1,0, 0, 1)1 : 100% right circular polarization (i.e., negative helicity); 

• (1, 0, 0, —1)/ : 100% left circular polarization (i.e., positive helicity). 

24.3 Relation Between Stokes Parameters of Incident and Scattered Radiation: 
The Mueller Matrix 

It is convenient to describe the scattering properties of a finite target in terms of the 4x4 Mueller ma- 
trix Sij relating the Stokes parameters (7^, Qi,Ui,Vi) and (I S ,Q S , U S ,V S ) of the incident and scattered 
radiation: 
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Once the amplitude scattering matrix elements are obtained, the Mueller matrix elements can be com- 
puted (Bohren & Huffman 1983): 
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= (|Si| 2 + |S 2 
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= (\Si\ 2 + \S 2 


l 2 -l£ 3 l 2 - 


|S 4 | 2 ) /2 


S23 


= Re (S 2 S% - 


SiS* 4 ) , 
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= Im(S 2 SI + 


SiSl) , 
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= Rc(S2S* 4 + 
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= Im (S2SI + 
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= Im (S4S?, - 
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= hntfiSS- 
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= Rc^SZ- 


S 3 S 4 *) . 





These matrix elements are computed in DDSCAT and passed to subroutine WRITESCA which handles 
output of scattering properties. Although the Muller matrix has 16 elements, only 9 are independent. 



24 SCATTERING B Y FINITE TARGETS: THE MUELLER MATRIX 



64 



The user can select up to 9 distinct Muller matrix elements to be printed out in the output files 

wxxxryyykzzz . sea and vixxxryyyori . a vg; this choice is made by providing a list ofindicesin ddscat .par 
(see Appendix lAt. 

If the user does not provide a list of elements, WRITE SCA will provide a "default" set of 6 selected 
elements: Su, S21, S31, S41 (these 4 elements describe the intensity and polarization state for scattering 
of unpolarized incident radiation), Si 2, and S13. 

In addition, WRITE SCA writes out the linear polarization P of the scattered light for incident unpo- 
larized light (see Bohren & Huffman 1983): 

(c2 , o2 \l/2 

P= ^l+^l) _ (82) 

Su 

24.4 Polarization Properties of the Scattered Radiation 

The scattered radiation is fully characterized by its Stokes vector (I S ,Q S ,U S ,V S ). As discussed in 
Bohren & Huffman (eq. 2.87), one can determine the linear polarization of the Stokes vector by op- 
erating on it by the Mueller matrix of an ideal linear polarizer: 



e _ 1 



/ 1 cos2£ sin 2^ \ 

cos 2^ cos 2 2^ cos 2£ sin 2£ 

sin 2^ sin 2£ cos 2£ sin 2 2£ 

\ / 



(83) 



where £ is the angle between the unit vector 6 S parallel to the scattering plane ("SP") and the "transmis- 
sion" axis of the linear polarizer. Therefore the intensity of light polarized parallel to the scattering plane 
is obtained by taking £ = and operating on (I s , Q S ,U S ,V S ) to obtain 

7(E S ||SP) = ^(h + Qs) (84) 

= [(Su + S 21 )h + (S 12 + S 22 )Q l + (S 13 + S 23 )Ui + (Su + S 2i )Vi] . (85) 

Similarly, the intensity of light polarized perpendicular to the scattering plane is obtained by taking 

£ = tt/2: 

/(E S ±SP) = \(I a -Qs) (86) 
[{S n - S 21 )Ii + (S12 - S 22 )Q 4 + (Si 3 - S 23 )C/, + (S14 - S 24 )y 4 ] . (87) 



2yfc 2 r 2 

24.5 Relation Between Mueller Matrix and Scattering Cross Sections 

Differential scattering cross sections can be obtained directly from the Mueller matrix elements by noting 
that 

l_ ( dC sl 
r 2 { dtl 



I s = — \ -17T ) h ■ (88) 



S.l 



Let SP be the "scattering plane": the plane containing the incident and scattered directions of propaga- 
tion. Here we consider some special cases: 

• Incident light unpolarized: Stokes vector s, = 1(1, 0, 0, 0): 

- cross section for scattering with polarization E s || SP: 

^ - 5P (!»!• + 1*!")- 51 <a. + ao m 

- cross section for scattering with polarization E s _L SP: 

iff -5^ (1*1* + !&!•)- <m 
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- total intensity of scattered light: 

= ^ W + 1*1' + l^l 2 + l^l 2 ) = ±S U (91) 

• Incident light polarized with Ej || SP: Stokes vector Sj = J(l, 1, 0, 0): 

- cross section for scattering with polarization E s || to SP: 

— = ^|S2| 2 = ^(Sii+Si2 + S2i+S2 2 ) (92) 



dft /t 2 ' 1 2fc 
- cross section for scattering with polarization E s _L to SP: 



L '<? 4 | 2 = -i ? (5 11 + 5 12 -5 21 -5 22 ) (93) 



dSl fc 2 ' 1 2fc 2 
total scattering cross section: 



a =^(|5 2 | 2 + |5 4 | 2 ) = ^(5 11 + 5 12 ) (94) 



dn k 2 Vl 1 1 1 y fc 2 

• Incident light polarized with Ej _L SP: Stokes vector Sj = 7(1, —1, 0, 0): 
- cross section for scattering with polarization E s || to SP: 



dC sca _ J_i 1 2 _ _j_ 
dSl ~ #»' 81 ~ 2fc 2 



72 l S 3| 2 = 7T72 (Sn -S 12 + S 2 i -fe) (95) 



- cross section for scattering with polarization E s _L SP: 

dC sca 1 i si i9 1 



- jfe2i fif i| 2 = 2^ (Sii - Sia - S21 + S22) (96) 

- total scattering cross section: 

^ = ^(\S^ + \S 3 f) = ^( Sll -S 12 ) (97) 

• Incident light linearly polarized at angle 7 to SP: Stokes vector Si = 7(1, cos 27, sin 27, 0): 

- cross section for scattering with polarization E s || to SP: 

dC 1 

~dW = 2F [(5n + 52l) + ( ^ 12 + 522) cos27 + + s ' 23) sill27] (98) 

- cross section for scattering with polarization E s _L SP: 

= i [{Sn ~ s ' 2i) + ( ^ 12 ~ S22) cos27 + ( ^ 13 - s ' 23) sin27] (99) 

- total scattering cross section: 

dC sca _ 1 



7T [5 , ii+5i 2 cos2 7 + 5i 3 sm27] (100) 
as I k z 

24.6 One Incident Polarization State Only (1 0RTH= 1) 

In some cases it may be desirable to limit the calculations to a single incident polarization state - for 
example, when each solution is very time-consuming, and the target is known to have some symmetry 
so that solving for a single incident polarization state may be sufficient for the required purpose. In this 
case, set IORTH=l in ddscat .par. 

When I0RTH=1, only fu and / 21 are available; hence, DDSCAT cannot automatically generate 
the Mueller matrix elements. In this case, the output routine WRITE SCA writes out the quantities |/n| 2 , 
|/ 2 i| 2 , Re (/1 1/21), and Im(/n/ 21 ) for each of the scattering directions. 
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The differential scattering cross section for scattering with polarization E s || and _L to the scattering 
plane are 

(§f),,, - f""' 2 < 101 » 
(tstL = ^ hlf <102) 

Note, however, that if IPHI is greater than 1, DDSCAT will automatically set I0RTH=2 even if 
ddscat . par specified I0RTH=1: this is because when more than one value of the target orientation 
angle $ is required, there is no additional "cost" to solve the scattering problem for the second incident 
polarization state, since when solutions are available for two orthogonal states for some particular target 
orientation, the solution may be obtained for another target orientation differing only in the value of 
$ by appropriate linear combinations of these solutions. Hence we may as well solve the "complete" 
scattering problem so that we can compute the complete Mueller matrix. 



25 Scattering by Periodic Targets: Generalized Mueller Matrix 

The Mueller scattering matrix Sij(9) described above was originally defined (see, e.g., Bohren & Huff- 
man 1983) to describe scattering of incident plane waves by finite targets, such as aerosol particles or 
dust grains. 

Draine & Flatau (2008) have extended the Mueller scattering matrix formalism to also apply to 
targets that are periodic and infinite in one or two dimensions (e.g., an infinite chain of particles, or a 
two-dimensional array of particles). 

25.1 Mueller Matrix S^j for Targets Periodic in One Direction 

The dimensionless 4x4 matrix (M, Q describes the scattering properties of targets that are periodic 
in one dimension. In the radiation zone, for incident Stokes vector L m , the scattered Stokes vector in 
direction (M, Q (see 321 .21 is (see Draine & Flatau 2008, eq. 63) 

I sc ,,i(M,0 = YRzZ S S d) ( M 'OIin, j , (103) 

° 3=1 

where R is the distance from the target repetition axis. 

DDSCAT 7.1 reports the scattering matrix elements in the output files waaarbbbkccc . sea. 

25.2 Mueller Matrix s\f ] (M, N) for Targets Periodic in Two Directions 

For targets that are periodic in two directions, the scattering intensities are described by the dimensionless 
4x4 matrix Sj^ (M, iV). In the radiation zone, for incident Stokes vector I m , the scattered Stokes 
vector in scattering order (M, N) is (see Draine & Flatau 2008, eq. 64) 

4 

I sca>i (M,N) = J2S% d) (M,N)I ill , j , (104) 

3=1 

where integers (M, N) define the scattering order (see 321.31 ). For given (M, N) there are actually two 
possible scattering directions - corresponding to transmission and reflection (see 321.31 ). and for each 
there is a value of S^ d) (M, N). 

For targets with 2-d periodicity, the scattering matrix elements S^ d '(M, N) are directly related to 
the usual transmission coefficient and reflection coefficient (see Draine & Flatau 2008, eq. 69-71). 
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26 Composite Targets with Anisotropic Constituents 

Section[l9]includes targets composed of anisotropic materials (ANIELLIPS, ANIRCTNGL, ANI_ELL_2, 
ANI_ELL_3). However, in each of these cases it is assumed that the dielectric tensor of the target ma- 
terial is diagonal in the "Target Frame". For targets consisting of a single material (in a single domain), 
it is obviously possible to choose the "Target Frame" to coincide with a frame in which the dielectric 
tensor is diagonal. 

However, for inhomogeneous targets containing anisotropic materials as, e.g., inclusions, the optical 
axes of the constituent material may be oriented differently in different parts of the target. In this case, it 
is obviously not possible to choose a single reference frame such that the dielectric tensor is diagonalized 
for all of the target material. 

To extend DDSCAT to cover this case, we must allow for off-diagonal elements of the dielectric 
tensor, and therefore of the dipole polarizabilities. It will be assumed that the dielectric tensor e is 
symmetric (this excludes magnetooptical materials - check this). 

Let the material at a given location in the grain have a dielectric tensor with complex eigenvalues 
j = 1 — 3. These are the diagonal elements of the dielectric tensor in a frame where it is diagonalized 
(i.e., the frame coinciding with the principal axes of the dielectric tensor). Let this frame in which the 
dielectric tensor is diagonalized have unit vectors corresponding to the principal axes of the dielectric 
tensor. We need to describe the orientation of the "Dielectric Frame" (DF) - defined by the "dielectric 
axes" ej - relative to the Target Frame. It is convenient to do so with rotation angles analogous to the 
rotation angles 0, <i>, and j3 used to describe the orientation of the Target Frame in the Lab Frame: 
Suppose that we start with unit vectors e.j aligned with the target frame ij . 

1. Rotate the DF through an angle 6r>F around axis ytf, so that 8bf is now the angle between ei 
and xtf ■ 

2. Now rotate the DF through an angle 0df around axis xtf, m such a way that e 2 remains in the 
xtf — ei plane. 

3. Finally, rotate the DF through an angle /?df around axis ei. 

The unit vectors are related to the TF basis vectors xtf, YtFj %tf by ; 

ei = xtf cos 6*df + Ytf sin 6*df cos 0df + ztf sin 6*df sin ^>df (105) 
e 2 = — xtf sin^DF cos/?df + yTF[cos 6*df cos/3df cos(/)df — sin/3oF sin^DF] 

+ztf [cos 0df cos /3df sin c/)df + sin/3DF cos^df] (106) 
63 = xtf sin #df sin /3d f — Ytf [cos 6*df sin /3df cos </>df + cos /3df sin <^df] 

— ztf[cos 0df sin/3nF sin 0df — cos/3df cos^df] (107) 

or, equivalently: 

xtf = ei cos ^df ^ 62 sin ^df cos/3df + 63 sin^DF sin/3nF (108) 
yTF = ei sin 6*df cos c/>df + e 2 [cos #df cos /3df cos 4>r>F — sin /?df sin 4>~df] 

— e 3 [cos 6*df sin /3df cos 0df + cos /3df sin </>df] (109) 
ztf = ei sin #df sin </>df + e 2 [cos 6*df cos /3df sin </>df + sin /3df cos 4>df] 

—63 [cos 6>df sin /3df sin </>df — cos /?df cos </>df] (1 10) 

Define the rotation matrix = iq • : 

R 13 = (111) 

/ cos #df sin 8uf cos </>df sin 9of sin </>df \ 

I— sin#DF cos /3d f cos 9uf cos /3df cos 0df — sin/?DF sin 0df cos 9uf cos /3df sin <^DF+sin /3df cos <j>riF ] 
\ sin 9uf sin /3df — cos Sdf sin/3DF cos <^df — cos /3df sin 0df — cos 9df sin /3df sin 0df + cos /3df cos 0df/ 

(112) 

and its inverse 

(R-'h = (113) 
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(cos #df — sin #df cos /3df sin 6bf sin /3df \ 

sin 9df cos </>df cos #df cos /?df cos (/>df — sin /3d f sin ^>df — cos #df sin /3df cos 0df — cos /3d f sin 0df I 
sin Oof sin 0df cos #df cos /3df sin 0df + sin /3df cos 0df — cos #df sin /3df sin 0df + cos /3d f cos </>df/ 

(114) 

The dielectric tensor is diagonal in the DF. If we calculate the dipole polarizabilility tensor in the DF it 
will also be diagonal, with elements 



(115) 

The polarizability tensor in the TF is given by 

(a TF )„„ = RijioP^jkiR-^km (116) 




Thus, we can describe a general anisotropic material if we provide, for each lattice site, the three 
diagonal elements e? F of the dielectric tensor in the DF, and the three rotation angles #df, /?df, and 
0df- We first use the LDR prescription and the elements to obtain the three diagonal elements a£ F of 
the polarizability tensor in the DF. We then calculate the polarizability tensor a TF using eq. dl 16b . 



27 Postprocessing: Calculating E and B In or Near the Target 

For some scientific applications (e.g., surface-enhanced Raman scattering) one wishes to calculate the 
electric field just outside the target boundary. It may also be of interest to calculate the electromagnetic 
field at positions within the target volume. DDSCAT includes tools to carry out such calculations. 

When DDSCAT is run with option IWRP0L=1, it creates files vrxxxryyykzzz . poln for «=1 and (if 
I0RTH=2) n=2. These files store the solution for the dipole polarizations, allowing postprocessing. Of 
particular interest is calculation of the electromagnetic field at locations in or near the target. 

To facilitate this, we provide a Fortran subroutine READPOL (in the file readpol . f 90) that reads 
in data from the wxxxryyy~kzzz .poln file. 

In addition, we provide a Fortran program DDfield.f90 that calls subroutine READPOL, and uses 
the data read in by READPOL to calculate the electric and magnetic field at locations specified by the 
user. The E and B contributions from the oscillating dipoles are calculated [including retardation effects 
- see, e.g., equations (7) and (10) from Draine & Weingartner (1996)] giving the "exact" electromagnetic 
field contributed by the oscillating dipoles. The incident plane wave field is added to the field generated 
by the oscillating dipoles to obtain the total E and B at the specified location. 

The user is interested in the continuous and finite field in a real target, as opposed to the field gener- 
ated by an array of oscillating point dipoles, which diverges as one approaches any one of the dipoles. 
To suppress the divergences, the contributions of dipoles within a distance d of the selected point is ob- 
tained by multiplying the "exact" contribution from the point dipole by a factor </>(r) that goes to zero 
more rapidly than r 3 , thus ensuring that the electric field contribution from an individual dipole goes to 
zero as we approach the dipole location. This is desired, because the discrete dipole approximation is 
based on the assumption that a given dipole is polarized by the electric field from the incident wave plus 
the electric field contributed by all other dipoles. When dealing with periodic (i.e., infinite) targets, the 
function </>(r) smoothly suppresses the contributions from distant dipoles, allowing the summations to be 
truncated. Please refer to Draine & Flatau (2008) for further discussion of (f>(r) and examples of results 
calculated by D D f i e 1 d. 

To calculate the electric and magnetic field at user-selected locations: 

• Run ddscat with WRPOL=l to create a file with the stored polarization (for example, 

wOOOrOOOkOOO .poll) 

• Compile and link the DD field code, by entering the src directory and typing 

make ddfield 

This will create an executable named DDf ield 
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• Use an editor to create a file ddfield. in providing the name of the file (in quotes!, where 
"quote " should be the "apostrophe " character) containing the stored target polarization, 
followed by the value of the interaction cutoff parameter 7 (see Draine & Flatau 2008, eq. 79) to 
be used, 

followed by the information defining a "track" along which E and B will be sampled. The track 

is specified by the beginning and end points (x±, y±, z\) and (xjv, Un, zn) in the TF, in units of 

dipole spacing d, followed by an integer specifying the number N of points in the track. Thus, for 

example: 

'filename' 

7 

xi/d yi/d z\/d Xfj/d vn/A z^/d N for track 1 
X\/d yi/d z\/d xn/cI vn/A z^/d N for track 2 

x\jd yi/d zi/d XN/d yw/d zjv/rf N for last track 

Here Xj/d yj/d zj/d are integer or decimal numbers giving the coordinates (in the "Tar- 
get Frame", in units of the lattice spacing d) of the positions at which the electromagnetic field is 
to be evaluated. 

• Type 

ddfield 
to run the executable. 

• The output will be written to two files, DDfield.E and DDBfield.B, containing the electric 
and magnetic fields (given as complex vectors) at each chosen location. The fields are evaluated at 
a time when the factor e~ tuJt is 1 (e.g., t — 0, 2ir/ui, ....). 

27.1 Sample Calculation 

Here is an example of an input file ddfield. in: 

' wOO Or 00 OkOO . poll ' = name of file with stored polarization 

5.00e-3 = gamma (interaction cutoff parameter) 

-64. 0. 0. 32. 0. 0. 481 (xa,ya,za , xb,yb,zb , NAB) 

The above example will sample the EM field at 48 1 points uniformly spaced along a track running 
from (x/d, y/d,z/d) = (-64, 0, 0) to (32, 0, 0) at intervals of Sx = 0.2d. 
Here's a second example of an input file ddfield . in: 

' wOO Or 00 OkOO . poll ' = name of file with stored polarization 

5.00e-3 = gamma (interaction cutoff parameter) 

-64. 32.64 32.64 32 32.64 32.64 481 (xa,ya,za , xb,yb,zb , NAB) 

Note that x, y, z are specified in the Target Frame, and in units of d. To properly locate these you 
will need to know the location of the TF origin (0, 0, 0) - see the target descriptions in ^191201 

We have run DDf ield on the output created by the "test problem" of a rectangular gold block with 
dimensions 0.5/im x 1/zm x ljum illuminated by a plane wave with wavelength A = O.bfim incident 
normal to the x 1/zm surface. The Au target has refractive index m — 0.9656 + 1.8628i (dielectric 
function e = -2.5374 + 3.5975i). (see ®. 

For the {x, 0, 0) track running through the center of the slab, the output DDf ield . E file will look 

like 

131072 = number of dipole s in Target 
Extent of occupied lattice sites 
1 32 = JXMIN, JXMAX 

1 64 = JYMIN, JYMAX 

1 64 = JZMIN, JZMAX 

-32.000000 0.000000 = (x_TF/d) min, (x_TF/d) max 

-32.000000 32.000000 = ( y_TF/d) min, ( y_TF/d) max 
-32.000000 32.000000 = ( z_TF/d) min , ( z_TF/d) max 
-0.499997 0.000000 = xmin ( TF ) , xmax ( TF ) (phys . units) 

-0.499997 0.499997 = ymin ( TF ) , ymax ( TF ) (phys. units) 

-0.499997 0.499997 = zmin ( TF ) , zmax ( TF ) (phys. units) 

0.015625 = d (phys units) 
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Fi gure 1 1 ! Solid line: \E\ 2 on track passing through center of Au block for the sample problem calculated using target option 
RCTGLPRSM. Dotted lines show where the Au block surface is located. Note how weak the E field is inside the Au. Dashed line: 
\E\ 2 on track passing near the corner of the slab. Incident radiation is propagating in the +:r direction. 



0.000000 0.000000 = PYD,PZD = period_y/dy, period_z/dz 
0. 000000 0.000000 = period_y, period_z (phys . units) 

0.500000 = wavelength in ambient medium (phys units ) 

0.196349 = k_x*d for incident wave 

0.000000 = k_y*d for incident wave 

0.000000 = k_z*d for incident wave 



5 . OOOE-03 = gamma (parameter for summation cutoff) 
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In the example given, the incident wave is propagating in the +x direction. The first layer of dipoles 
in the rectangular target is at x/d = —31.5, and the last layer is at x/d = —0.5. The "surface" of the 
target is at x/d = —32 (x = — 0.50/im) and x/d = (x = 0). 

Figure fTTIshows how \E\ 2 varies along the tracks specified by the above two sample ddf ield . in 
files. As expected, the wave is strongly suppressed in the interior of the Au block. 
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28 Finale 

This User Guide is somewhat inelegant, but we hope that it will prove useful. The structure of the 
ddscat .par file is intended to be simple and suggestive so that, after reading the above notes once, 
the user may not have to refer to them again. 

Known bugs in DDSCAT will be posted at the DDSCAT web site, 

|http : / / www .astro . princeton . edu /( ^draine /DDSCAT . html 
and the latest version of DDSCAT will always be available at this site. The source code is also available 
at 

fFittp : //code . google . com/p/ddscat/| 
There is also a wiki: 

jhttp : //ddscat . wikidot . com/start| 

Users are encouraged to provide B. T. Draine (draine@astro.princeton.edu) with their 
email address; email notification of bug fixes, and any new releases of DDSCAT, will be made known 
to those who do! 

P. J. Flatau maintains the WWW page "SCATTERLIB - Light Scattering Codes Library" with URL 
|http ://atol.ucsd. edu/^pf latau The SCATTERLIB Internet site is a library of light scat- 
tering codes. Emphasis is on providing source codes (mostly FORTRAN). However, other information 
related to scattering on spherical and non-spherical particles is collected: an extensive list of references 
to light scattering methods, refractive index, etc. This URL page contains section on the discrete dipole 
approximation. 

Concrete suggestions for improving DDSCAT (and this User Guide) are welcomed. 

Users of DDSCAT may wish to cite papers describing DDA theory and its implementation, such as 
Draine (1988), Draine & Flatau (2004), and Draine & Flatau (2008). The latter reference describes ex- 
tension of the DDA to periodic structures, and generalization of the Mueller scattering matrix to describe 
scattering by 1-d and 2-d periodic structures. 

To cite this User Guide, we suggest the following citation: 

Draine, B.T., & Flatau, RJ. 2010, "User Guide for the Discrete Dipole Approximation Code 
DDSCAT (Version 7 . 1 )", |http://arxiv.org/abs/l 002. 1505| 

Finally, the authors have one special request: We would appreciate preprints and (especially ! ) reprints 
of any papers which make use of DDSCAT! 
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A Understanding and Modifying dd scat. par 

In order to use DDSCAT to perform the specific calculations of interest to you, it will be necessary to 
modify the ddscat . par file. Here we list the sample ddscat . par file, followed by a discussion of 
how to modify this file as needed. Note that all numerical input data in DDSCAT is read with free-format 
READ (IDEV, *) . . . statements. Therefore you do not need to worry about the precise format in which 
integer or floating point numbers are entered on a line. The crucial thing is that lines in ddscat . par 
containing numerical data have the correct number of data entries, with any informational comments 
appearing after the numerical data on a given line. 

' ========= Parameter file for v7 . 1 . ===================' 

'**** Preliminaries ****' 

' NOTORQ' = CMT0RQ*6 (NOTORQ, DOTORQ) — either do or skip torque calculations 

' PBCGS2 ' = CMDS0L*6 (PBCGS2 , PBCGST, PETRKP) — select solution method 

' GPFAFT ' = CMDFFT* 6 (GPFAFT, FFTMKL) 

' GKDLDR' = CALPHA*6 (GKDLDR, LATTDR 

'NOTBIN' = CBINFLAG (NOTBIN, ORIBIN, ALLBIN) 

'**** Initial Memory Allocation ****' 

100 100 100 = dimensioning allowance for target generation 
'**** Target Geometry and Composition ****' 
' RCTGLPRSM' = CSHAPE*9 shape directive 
32 64 64 = shape parameters 1-3 

1 = NCOMP = number of dielectric materials 
' . . /diel/Au_evap' = file with refractive index 1 
'**** Error Tolerance ****' 

1.00e-5 = TOL = MAX ALLOWED (NORM OF | G>=AC I E>-ACA | X> ) / (NORM OF AC|E>) 

'**** Interaction cutoff parameter for PBC calculations ****' 

1.00e-2 = GAMMA (le-2 is normal, 3e-3 for greater accuracy) 

'**** Angular resolution for calculation of <cos>, etc. ****' 

0.5 = ETASCA (number of angles is proportional to [ ( 3+x) /ETASCA] * 2 ) 

'**** Wavelengths (micron) ****' 

0.5000 0.5000 1 'LIN' = wavelengths ( first , last , how many , how=LIN, INV, LOG) 
'**** Effective Radii (micron) **** ' 

0.49237 0.49237 1 'LIN' = aeff ( first , last , how many , how=LIN, INV, LOG) 
'**** Define Incident Polarizations ****' 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = IORTH (=1 to do only pol. state eOl; =2 to also do orth. pol . state) 
'**** Specify which output files to write ****' 

1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient. 

1 = IWRPOL (=0 to suppress, =1 to write ".pol" file for each (BETA, THETA) 
'**** Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 

0. 0. 1 = THETMI, THETMX, NTHETA (theta=angle between al and k) 

0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 

'**** Specify first IWAV, IRAD, IORI (normally 0) ****' 

= first IWAV, first IRAD, first IORI (0 to begin fresh) 

'**** Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 
11 12 21 22 31 41 = indices ij of elements to print 
'**** Specify Scattered Directions ****' 

' LFRAME ' = CMDFRM (LFRAME, T FRAME for Lab Frame or Target Frame) 

2 = number of scattering planes 

0. 0. 180. 5 = phi, thetan_min, thetan_max, dtheta (in deg) for plane 1 
90. 0. 180. 5 = phi, thetan_min, thetan_max, dtheta (in deg) for plane 2 
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Lines Comments 
1-2 comment lines 

3 NOTORQ if torque calculation is not required; 
DOTORQ if torque calculation is required. 

4 PBCGS2 is recommended; other options are PBCGST and PETRKP (see ^TTb . 

5 GPFAFT is supplied as default, but FFTMKL is recommended if DDSCAT has been compiled with 
the Intel® Math Kernal Library (see 53631 [Hi 

6 GKDLDR is recommended as the prescription for polarizabilities (see Sjl3]) 

7 NOTB IN for no unformatted binary output. 

ORIBIN for unformatted binary dump of orientational averages only; 
ALLBIN for full unformatted binary dump ( £110.21 ); 

8 comment line 

9 initial memory allocation NX,NY,NZ. These must be large enough to accomodate the target that 
will be generated. 

10 comment line 

1 1 specify choice of target shape (see Sjl9]for description of options RCTGLPRSM, ELLIPSOID, 
TETRAHDRN, ...) 

12 shape parameters SHPAR1, SHPAR2, SHPAR3, ... (see g9]l. 

13 number of different dielectric constant tables (CHI). 

14 name(s) of dielectric constant table(s) (one per line). 

15 comment line 

16 TOL = error tolerance h: maximum allowed value of \A^E - A^AP\/\A^E\ [see eq.lQB]. 

17 comment line 

18 GAMMA = interaction cutoff parameter 7 (see Draine & Flatau 2009) 
the value of 7 does not affect calculations for isolated targets 

19 comment line 

20 ETASCA - parameter 77 controlling angular averages (<f23]). 

2 1 comment line 

22 A - first, last, how many, how chosen. 

23 comment line 

24 a c ff - first, last, how many, how chosen. 

25 comment line 

26 specify x,y,z components of (complex) incident polarization §01 ( ^22l) 

27 IORTH = 2 to do both polarization states (normal); 
IORTH = 1 to do only one incident polarization. 

28 comment line 

29 IWRKSC = to suppress writing of ".sea" files; 
IWRKSC = 1 to enable writing of ".sea" files. 

30 IWRPOL = to suppress writing of ".pol" files. 
IWRPOL = 1 to write ".pol" files. 

3 1 comment line 

32 j3 (see jfTTb - first, last, how many . 

33 (see £1171 - first, last, how many. 

34 $ (see qlTt - first, last, how many. 

35 comment line 

36 IWAVO IRAD0 IORI0 - starting values of integers IWAV IRAD I ORI (normally 0). 

37 comment line 

38 Ns = number of scattering matrix elements (must be < 9) 

39 indices ij of N$ elements of the scattering matrix Sly 

40 comment line 

41 specify whether scattered directions are to be specified by CMDFRM= ' L FRAME ' or ' T FRAME ' . 

42 = number of scattering planes to follow 

43 4> s for first scattering plane, 6 S)m i n , s>max , how many 6 S values; 
44,... (f> s for 2nd,... scattering plane, ... 
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B vrxxxryyy .avg Files 

The file wOOOrOOOori. avg contains the results for the first wavelength (wOOO) and first target radius 
(rOOO) averaged over orientations ( . avg). The wOOOrOOOori . avg file generated by the sample 
calculation should look like the following: 

DDSCAT DDSCAT 7.1.0 [10.02.04] 

TARGET Rectangular prism; NX,NY,NZ= 32 64 64 

GPFAFT method of solution 

GKDLDR prescription for polarizabilies 

RCTGLPRSM shape 

131072 = NATO = number of dipoles 

0.03173411 = d/aeff for this target [d=dipole spacing] 

0.015625 = d (physical units) 
AEFF= 0.492370 = effective radius (physical units) 

WAVE= 0.500000 = wavelength (physical units) 

K*AEFF= 6.187304 = 2*pi*aeff /lambda 

n= ( 0.9656 , 1.8628), eps.= ( -2.5374 , 3.5975) |m|kd= 0.4120 for subs. 1 

TOL= 1.000E-05 = error tolerance for CCG method 
( 1.00000 0.00000 0.00000) = target axis Al in Target Frame 
( 0.00000 1.00000 0.00000) = target axis A2 in Target Frame 

NAVG= 2292 = (theta,phi) values used in comp. of Qsca,g 
( 0.19635 0.00000 0.00000) = k vector (latt. units) in Lab Frame 
( 0.00000, 0.00000) ( 1.00000, 0.00000) ( 0.00000, 0. 00000) =inc.pol.vec. 1 in LF 
( 0.00000, 0.00000) ( 0.00000, 0.00000) ( 1.00000, 0.00000)=inc.pol.vec. 2 in LF 

0.000 0.000 = beta_min, beta_max ; NBETA = 1 

0.000 0.000 = theta_min, theta_max; NTHETA- 1 

0.000 0.000 = phi_min, phi_max ; NPHI = 1 



0.5000 = ETASCA = param. controlling # of scatt . dirs used to calculate <cos> etc. 
Results averaged over 1 target orientations 

and 2 incident polarizations 
Qext Qabs Qsca g(l)= ; <cos> <cos"2> Qbk Qpha 

JO=l: 3.6494E+00 1.2208E+00 2.4287E+00 4.6029E-01 8.7798E-01 2.6214E+00 -1.7030E-01 
JO=2: 3.6494E+00 1.2208E+00 2.4287E+00 4.6029E-01 8.7798E-01 2.6214E+00 -1.7030E-01 
mean: 3.6494E+00 1.2208E+00 2.4287E+00 4.6029E-01 8.7798E-01 2.6214E+00 -1.7030E-01 
Qpol= -4.5300E-06 dQpha= -1.0133E-06 

Qsca*g(l) Qsca*g(2) Qsca*g(3) iter mxiter Nsca 
JO=l: 1.1179E+00 -9.4339E-08 3.3191E-07 27 117 2292 

JO=2: 1.1179E+00 1.9882E-07 -2.6904E-08 29 117 2292 

mean: 1.1179E+00 5.2242E-08 1.5250E-07 

Mueller matrix elements for selected scattering directions in Lab Frame 
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223E- 


05 


-1 


148E- 


04 


40 


00 





00 





69023 


3 


. 7585E+01 


-2 


5942E+01 


-2 


594E+01 


3 


759E+01 


-5 


071E- 


05 


-9 


225E- 


05 


45 


00 





00 





93616 


2 


. 2608E+01 


-2 


1165E+01 


-2 


116E+01 


2 


261E+01 


-4 


507E- 


05 


-2 


219E- 


07 


50 


00 





00 





55260 


1 


. 6124E+01 


-8 


9100E+00 


-8 


910E+00 


1 


612E+01 


-2 


009E- 


05 


5 


881E- 


05 


55 


00 





00 





04414 


1 


. 5380E+01 


-6 


7895E-01 


-6 


789E-01 


1 


538E+01 


3 


102E- 


05 


1. 


028E- 


04 


60 


00 





00 





00152 


1 


. 4134E+01 


2 


1518E-02 


2 


152E-02 


1 


413E+01 


4 


444E- 


05 


9 


157E- 


05 


65 


00 





00 





25433 


1 


. 0795E+01 


-2 


7454E+00 


-2 


745E+00 


1 


079E+01 


5 


736E- 


05 


6 


703E- 


05 


70 


00 





00 





61640 


7 


. 1817E+00 


-4 


4268E+00 


-4 


427E+00 


7 


182E+00 


4 


303E- 


05 


3 


109E- 


05 


75 


00 





00 





84301 


4 


. 7205E+00 


-3 


9795E+00 


-3 


979E+00 


4 


721E+00 


2 


952E- 


05 


5 


581E- 


06 


80 


00 





00 





83319 


3 


. 3467E+00 


-2 


7884E+00 


-2 


788E+00 


3 


347E+00 


1 


446E- 


05 


-5 


168E- 


06 


85 


00 





00 





83162 


2 


. 5365E+00 


-2 


1094E+00 


-2 


109E+00 


2 


536E+00 


1 


247E- 


05 


-1 


128E- 


05 


90 


00 





00 





95712 


2 


. 1245E+00 


-2 


0334E+00 


-2 


033E+00 


2 


124E+00 


1 


498E- 


05 


-9 


050E- 


06 


95 


00 





00 





97305 


2 


. 1420E+00 


-2 


0843E+00 


-2 


084E+00 


2 


142E+00 


2 


119E- 


05 


-6 


176E- 


06 


100 


00 





00 





86407 


2 


. 4048E+00 


-2 


0779E+00 


-2 


078E+00 


2 


405E+00 


2 


250E- 


05 


-4 


506E- 


07 


105 


00 





00 





86090 


2 


. 6074E+00 


-2 


2447E+00 


-2 


245E+00 


2 


607E+00 


1 


945E- 


05 


6 


889E- 


07 


110 


00 





00 





95634 


2 


. 8343E+00 


-2 


7105E+00 


-2 


711E+00 


2 


834E+00 


1 


241E- 


05 


-5 


764E- 


06 


115 


00 





00 





84227 


3 


. 7313E+00 


-3 


1428E+00 


-3 


143E+00 


3 


731E+00 


7 


830E- 


06 


-2 


114E- 


05 


120 


00 





00 





53216 


5 


. 8481E+00 


-3 


1121E+00 


-3 


112E+00 


5 


848E+00 


4 


552E- 


06 


-4 


294E- 


05 


125 


00 





00 





31008 


8 


. 7198E+00 


-2 


7039E+00 


-2 


704E+00 


8 


720E+00 


1 


571E- 


05 


-5 


892E- 


05 
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130 


00 


.00 





21846 


1 


0852E+01 


-2 


3708E+00 


-2 


371E+00 


1 


085E+01 


3 


464K- 


05 


-6 


386E- 


05 


135 


00 


.00 





19080 


1 


0801E+01 


-2 


0610E+00 


-2 


061E+00 


1 


080E+01 


4 


566E- 


05 


-4 


999E- 


05 


140 


00 


.00 





14140 


8 


1464E+00 


-1 


1519E+00 


-1 


152E+00 


8 


146E+00 


4 


584E- 


05 


-2 


501E- 


05 


145 


00 


.00 





01267 


3 


8742E+00 


4 


9092E-02 


4 


909E-02 


3 


874E+00 


3 


123E- 


05 


-1 


770E- 


06 


150 


00 


.00 





35726 


1 


7618E+00 


-6 


2942E-01 


-6 


294E-01 


1 


762E+00 


5 


787E- 


06 


1 


075E- 


05 


155 


00 


.00 





45976 


: 


1004E+01 


-5 


0593E+00 


-5 


059E+00 


: 


100E+01 


-2 


023E- 


05 


6 


036E- 


06 


160 


00 


.00 





25129 


4 


5243E+01 


-1 


1369E+01 


-1 


137E+01 


4 


524E+01 


-3 


670E- 


05 


-1 


867E- 


05 


165 


00 


.00 





12705 


i 


1280E+02 


-1 


4331E+01 


-1 


433E+01 


1 


128E+02 


-2 


745E- 


05 


-4 


470E- 


05 


170 


00 


.00 





05271 


2 


0294E+02 




0697E+01 


-1 


070E+01 


2 


029E+02 


-8 


684E- 


06 


-9 


197E- 


05 


175 


00 


.00 





01266 


2 


8300E+02 


-3 


5830E+00 


-3 


583E+00 


2 


830E+02 


6 


964E- 


06 


-1 


170E- 


04 


180 


00 


.00 





00000 


3 


1528E+02 


1 


0529E-03 


1 


053E-03 


3 


153E+02 


3 


305E- 


05 


-1 


139E- 


04 





00 


90 .00 





00000 


1 


2306E+03 


: 


4038E-03 


1 


404E-03 


: 


231E+03 


6 


016E- 


05 


-2 


801E- 


05 


5 


00 


90 .00 





01833 


1 


0648E+03 


-l 


9514E+01 


-1 


951E+01 


: 


065E+03 


2 


780E- 


05 


-3 


678E- 


05 


10 


00 


90 .00 





08033 


6 


7623E+02 


-5 


4322E+01 


-5 


432E+01 


6 


762E+02 


1 


102E- 


05 


-3 


231E- 


05 


15 


00 


90 .00 





21697 


2 


9299E+02 


-6 


3570E+01 


-6 


357E+01 


2 


930E+02 


1 


738E- 


05 


-2 


438E- 


05 


20 


00 


90 .00 





52318 


7 


4600E+01 


-3 


9029E+01 


-3 


903E+01 


7 


460E+01 


2 


034E- 


05 


-1 


601E- 


05 


25 


00 


90 .00 





47167 


2 


1356E+01 




0073E+01 


-1 


007E+01 


2 


136E+01 


-1 


107E- 


06 


-1 


376E- 


05 


30 


00 


90 .00 





09869 


3 


8576E+01 


-3 


8071E+00 


-3 


807E+00 


3 


858E+01 


-2 


730E- 


05 


-1 


969E- 


05 


35 


00 


90 .00 





33427 


4 


8965E+01 


-1 


6367E+01 


-1 


637E+01 


4 


896E+01 


-4 


634E- 


05 


-2 


751E- 


05 


40 


00 


90 .00 





69023 


3 


7585E+01 


-2 


5943E+01 


-2 


594E+01 


3 


759E+01 


-3 


474E- 


05 


-2 


855E- 


05 


45 


00 


90 .00 





93616 


2 


2608E+01 


-2 


1165E+01 


-2 


116E+01 


2 


261E+01 


-2 


394E- 


06 


-2 


566E- 


05 


50 


00 


90 .00 





55260 


1 


6124E+01 


-8 


9099E+00 


-8 


910E+00 


1 


612E+01 


3 


199E- 


05 


-2 


475E- 


05 


55 


00 


90 .00 





04413 


1 


5380E+01 


-6 


7869E-01 


-6 


787E-01 


1 


538E+01 


5 


216E- 


05 


-3 


188E- 


05 


60 


00 


90 .00 





00154 


1 


4133E+01 


2 


1708E-02 


2 


171E-02 


1 


413E+01 


4 


590E- 


05 


-3 


833E- 


05 


65 


00 


90 .00 





25432 


1 


0795E+01 


-2 


7453E+00 


-2 


745E+00 


1 


079E+01 


2 


070E- 


05 


-3 


406E- 


05 


70 


00 


90 .00 





61638 


7 


1818E+00 


-4 


4267E+00 


-4 


427E+00 


7 


182E+00 


-6 


848E- 


06 


-1 


910E- 


05 


75 


00 


90 .00 





84300 


4 


7205E+00 


-3 


9794E+00 


-3 


979E+00 


4 


721E+00 


-2 


576E- 


05 


9 


569E- 


07 


80 


00 


90 .00 





83319 


3 


3466E+00 


-2 


7883E+00 


-2 


788E+00 


3 


347E+00 


-2 


862E- 


05 


1 


663E- 


05 


85 


00 


90 .00 





83162 


2 


5364E+00 


-2 


1093E+00 


-2 


109E+00 


2 


536E+00 


-1 


826E- 


05 


2 


157E- 


05 


90 


00 


90 .00 





95713 


2 


1244E+00 


-2 


0333E+00 


-2 


033E+00 


2 


124E+00 


-7 


100E- 


08 


1 


742E- 


05 


95 


00 


90 .00 





97305 


2 


1420E+00 


-2 


0842E+00 


-2 


084E+00 


2 


142E+00 


1 


977E- 


05 


1 


030E- 


05 


100 


00 


90 .00 





86407 


2 


4047E+00 


-2 


0778E+00 


-2 


078E+00 


2 


405E+00 


3 


345E- 


05 


7 


257E- 


06 


105 


00 


90 .00 





86090 


2 


6073E+00 


-2 


2446E+00 


-2 


245E+00 


2 


607E+00 


3 


586E- 


05 


1 


214E- 


05 


110 


00 


90 .00 





95634 


2 


8342E+00 


-2 


7104E+00 


-2 


710E+00 


2 


834E+00 


2 


374E- 


05 


2 


154E- 


05 


115 


00 


90 .00 





84225 


3 


7312E+00 


-3 


1426E+00 


-3 


143E+00 


3 


731E+00 


6 


297E- 


07 


2 


754E- 


05 


120 


00 


90 .00 





53214 


5 


8481E+00 


-3 


1120E+00 


-3 


112E+00 


5 


848E+00 


-2 


385E- 


05 


2 


262E- 


05 


125 


00 


90 .00 





31007 


8 


7197E+00 


-2 


7037E+00 


-2 


704E+00 


8 


720E+00 


-3 


821E- 


05 


6 


028E- 


06 


130 


00 


90 .00 





21846 


1 


0852E+01 


-2 


3708E+00 


-2 


371E+00 


1 


085E+01 


-3 


448E- 


05 


-1 


354E- 


05 


135 


00 


90 .00 





19080 


1 


0801E+01 


-2 


0610E+00 


-2 


061E+00 


1 


080E+01 


-1 


698E- 


05 


-2 


486E- 


05 


140 


00 


90 .00 





14140 


8 


1462E+00 


-1 


1519E+00 


-1 


152E+00 


8 


146E+00 


1 


131E- 


06 


-2 


232E- 


05 


145 


00 


90 .00 





01267 


3 


8740E+00 


4 


9075E-02 


4 


908E-02 


3 


874E+00 


6 


794E- 


06 


-9 


904E- 


06 


150 


00 


90 .00 





35722 


1 


7617E+00 


-6 


2933E-01 


-6 


293E-01 


1 


762E+00 


-1 


560E- 


07 


1 


926E- 


06 


155 


00 


90 .00 





45974 


1 


1005E+01 


-5 


0592E+00 


-5 


059E+00 


1 


100E+01 


-8 


574E- 


06 


1 


684E- 


06 


160 


00 


90 .00 





25128 


4 


5243E+01 


-1 


1369E+01 


-1 


137E+01 


4 


524E+01 


-7 


389E- 


06 


-1 


695E- 


05 


165 


00 


90 .00 





12704 


: 


1280E+02 




4331E+01 


-1 


433E+01 


: 


128E+02 


1 


518E- 


06 


-5 


206E- 


05 


170 


00 


90 .00 





05271 


2 


0294E+02 


-l 


0698E+01 


-1 


070E+01 


2 


029E+02 


8 


093E- 


06 


-9 


152E- 


05 


175 


00 


90 .00 





01267 


2 


8300E+02 


-3 


5861E+00 


-3 


586E+00 


2 


830E+02 


-2 


762E- 


06 


-1 


182E- 


04 


180 


00 


90 .00 





00000 


3 


1528E+02 




0834E-03 


-1 


083E-03 


3 


153E+02 


-3 


244E- 


05 


-1 


142E- 


04 
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C vsxxxryyy^zzz .sea Files 

The wOOOrOOOkOOO . sea file contains the results for the first wavelength (wOOO), first target radius 
(rOOO), and first orientation (kOOO). The wOOOrOOOkOOO.sca file created by the sample calculation 
should look like the following: 



DDSCAT DDSCAT 7.1.0 [10.02.04] 

TARGET Rectangular prism; NX,NY,NZ= 32 64 64 

GPFAFT method of solution 

GKDLDR prescription for polarizabilies 

RCTGLPRSM shape 

131072 = NATO = number of dipoles 

0.03173411 = d/aeff for this target [d=dipole spacing] 
0.015625 = d (physical units) 

physical extent of target volume in Target Frame 

-0.499997 0.000000 = xmin,xmax (physical units) 

-0.499997 0.499997 = ymin,ymax (physical units) 

-0.499997 0.499997 = zmin, zmax (physical units) 

AEFF= 0.492370 = effective radius (physical units) 

WAVE= 0.500000 = wavelength (physical units) 

K*AEFF= 6.187304 = 2*pi*aeff /lambda 

n= ( 0.9656 , 1.8628), eps . = ( -2.5374 , 3.5975) I m | kd 
TOL= 1.000E-05 = error tolerance for CCG method 
1.00000 0.00000 0.00000) = target axis Al in Target Frame 
0.00000 1.00000 0.00000) = target axis A2 in Target Frame 
NAVG= 2292 = (theta,phi) values used in comp. of Qsca,g 



0.4120 for subs. 1 



.19635 
.00000 
0.00000 
1.00000 
.00000 
.19635 



0.00000 0.00000) 
. 00000) ( 1 .00000, 
. 00000) ( .00000, 
0.00000 0.00000) 
1.00000 0.00000) 
0.00000 0.00000) 



k vector (latt. units) in TF 
0.00000) ( 0.00000, . 00000) =inc .pol .vec . 
0.00000) ( 1.00000, . 00000) =inc .pol .vec . 2 in TF 
target axis Al in Lab Frame 
target axis A2 in Lab Frame 
k vector (latt. units) in Lab Frame 
0.00000, 0.00000) ( 1.00000, 0.00000) ( 0.00000, 0.00000)=inc.pol.vec. 
0.00000, 0.00000) ( 0.00000, 0.00000) ( 1.00000, 0.00000)=inc.pol.vec. 
BETA = 0.000 = rotation of target around Al 
THETA= 0.000 = angle between Al and k 

PHI = 0.000 = rotation of Al around k 
0.5000 = ETASCA = param. controlling # of scatt . dirs used to calculate <cos> etc. 



1 in TF 



1 in LF 

2 in LF 



JO=l 
JO=2 
mean 
Qpol 

JO=l 
JO=2 
mean 



Qext Qabs Qsca 

3.6494E+00 1.2208E+00 2.4287E+00 

3.6494E+00 1.2208E+00 2.4287E+00 

3.6494E+00 1.2208E+00 2.4287E+00 
-4 . 5300E-06 

Qsca*g(l) Qsca*g(2) Qsca*g(3) 

1.1179E+00 -9.4339E-08 3.3191E-07 
1.1179E+00 1.9882E-07 -2.6904E-08 
1.1179E+00 5.2242E-08 1.5250E-07 



g ( 1 ) =<cos> 
4 . 6029E-01 
4 . 6029E-01 
4 . 6029E-01 



<cos"2> Qbk 
S.7798E-01 2.6214E+00 
S.7798E-01 2.6214E+00 
5.7798E-01 2.6214E+00 
dQpha= 



Qpha 
-1 . 7030E 
-1 . 7030E 
-1 . 7030E 
-1 . 0133E 



01 
01 
01 
-06 



iter 
27 
29 



mxiter 
117 
117 



Nsca 
2292 
2292 



Mueller matrix elements for selected scattering directions in Lab Frame 



theta 


phi 




Pol . 




3_11 




3_12 




3_21 




3_22 




3_31 






3_41 







00 





00 





00000 


1 


. 2306E + 03 


-1 


. 4038E-03 


-1 


. 404E-03 


1 


. 231E+03 


-5 


. 918E- 


05 


-2 


. 785E- 


05 


5 


00 





00 





01832 


1 


. 0648E+03 


-1 


. 9511E+01 


-1 


. 951E+01 


1 


. 065E+03 


-9 


. 344E- 


05 


8 


. 054E- 


07 


10 


00 





00 





08034 


6 


. 7623E+02 


-5 


. 4325E+01 


-5 


. 433E+01 


6 


. 762E+02 


-9 


. 059E- 


05 


1 


. 788E- 


05 


15 


00 





00 





21696 


2 


. 9300E+02 


-6 


. 3568E+01 


-6 


. 357E+01 


2 


. 930E+02 


-4 


. 868E- 


05 


7 


. 888E- 


05 


20 


00 





00 





52316 


7 


. 4600E+01 


-3 


. 9028E+01 


-3 


. 903E+01 


7 


. 460E+01 


1 


. 372E- 


05 


4 


. 646E- 


05 


25 


00 





00 





47166 


2 


. 1356E+01 


-1 


. 0073E+01 


-1 


. 007E+01 


2 


. 136E+01 


3 


. 462E- 


05 


-1 


. 389E- 


05 


30 


00 





00 





09869 


3 


. 8576E+01 


-3 


. 8070E+00 


-3 


. 807E+00 


3 


. 858E+01 


1 


. 067E- 


05 


-8 


. 347E- 


05 


35 


00 





00 





33427 


4 


. 8965E+01 


-1 


. 6368E+01 


-1 


. 637E+01 


4 


. 897E+01 


-2 


.223E- 


05 


-1 


. 148E- 


04 


40 


00 





00 





69023 


3 


. 7585E+01 


-2 


. 5942E+01 


-2 


. 594E+01 


3 


. 759E+01 


-5 


. 071E- 


05 


-9 


.225E- 


05 


45 


00 





00 





93616 


2 


.2608E+01 


-2 


. 1165E+01 


-2 


. 116E+01 


2 


.261E+01 


-4 


. 507E- 


05 


-2 


.219E- 


07 


50 


00 





00 





55260 


1 


. 6124E+01 


-8 


. 9100E+00 


-8 


. 910E+00 


1 


. 612E+01 


-2 


.009E- 


05 


5 


. 881E- 


05 


55 


00 





00 





04414 


1 


. 5380E+01 


-6 


. 7895E-01 


-6 


. 789E-01 


1 


. 538E+01 


3 


. 102E- 


05 


1 


. 028E- 


04 


60 


00 





00 





00152 


1 


. 4134E+01 


2 


. 1518E-02 


2 


. 152E-02 


1 


. 413E+01 


4 


. 444E- 


05 


9 


. 157E- 


05 


65 


00 





00 





25433 


1 


. 0795E+01 


-2 


. 7454E+00 


-2 


. 745E+00 


1 


. 079E+01 


5 


. 736E- 


05 


6 


. 703E- 


05 


70 


00 





00 





61640 


7 


. 1817E+00 


-4 


. 4268E+00 


-4 


. 427E+00 


7 


. 182E+00 


4 


.303E- 


05 


3 


. 109E- 


05 


75 


00 





00 





84301 


4 


. 7205E+00 


-3 


. 9795E+00 


-3 


. 979E+00 


4 


. 721E+00 


2 


. 952E- 


05 


5 


.581E- 


06 


80 


00 





00 





83319 


3 


. 3467E+00 


-2 


. 7884E+00 


-2 
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D wxxxryyy Yzzz . p o In Files 

The wOOOrOOOkOOO.poll file contains the polarization solution for the first wavelength ( w ), first 
target radius (r 0), first orientation (k 0), and first incident polarization (poll). In order to limit the 
size of this file, it has been written as an unformatted file. This preserves full machine precision for the 
data, is quite compact, and can be read efficiently, but unfortunately the file is not fully portable because 
different computer architectures (e.g., Linux vs. MS Windows) have adopted different standards for 
storage of "unformatted" data. However, anticipating that many users will be computing within a single 
architecture, the distribution version of DDSCAT uses this format. 

Additional warning: even on a single architecture, users should be alert to the possibility that different 
compilers may follow different conventions for reading/writing unformatted files. 

The file contains the following information: 

• The location of each dipole in the target frame. 

• (k x , k y , k z )d, where k is the incident k vector. 

• (Eo x , Eq v> Eoz), the complex polarization vector of the incident wave. 

• a -1 e? 3 , the inverse of the symmmetric complex polarizability tensor for each of the dipoles in the 
target. 

• (P x ,P y , P z ), the complex polarization vector for each of the dipoles. 

The interested user should consult the routine writepol . f, or the to see how this information has 
been organized in the unformatted file. The user can, of course, simply use the routine readpol . f to 
read this file and load the information into memory. 

Note that a subroutine readpol . f 90 is already provided that reads the stored file, and the program 
DD field, f 90 (which uses readpol . f 90) is provided to calculate E and B at user-selected postions 
-see $27] 
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Qabs, 11 
Qext, 1 1 
Qsca.i 1 1 
Qpha> 11 

Sj - scattering amplitude matrix, 62 

Sij - 4x4 Mueller scattering matrix, 61, 63, 64 

Sj - scattering amplitude matrix, 62 

i> - target orientation angle, 20, 28 

$df, 42, 67 

9 - target orientation angle, 20, 28 
Gdf, 42, 67 

f3 - target orientation angle, 20, 28 

0df, 42, 67 

rj - parameter for choosing scattering angles, 61 
77 - parameter for selection of scattering angles, 60 
ai, a 2 target vectors, 20 
<f> s - scattering angle, 55, 59 
9 S - scattering angle, 55, 59 
a eS , 7, 19 

fij - amplitude scattering matrix, 61 
fij - scattering amplitude matrix, 62 
Qr - radiative torque efficiency vector, 27 
Q pr - radiative force efficiency vector, 27 

absorption efficiency factor Q a b s , 1 1 
averages over scattering angles, 58 

CALLTARGET, 43 
CALLTARGET and PBC, 43 
CMDFFT - specifying FFT method, 74 
CMDFRM - specifying scattering directions, 55, 74 
compiling and linking, 14 
conjugate gradient algorithm, 23 

DDfield, 20, 68 

ddscat.logDOO - output file, 22 

ddscat.out, 22 

ddscat.par 

ALLBIN, 23 

CSHAPE, 32 

diel.tab, 18 

DOTORQ, 27 

ETASCA, 27 

GAMMA, 19 

GKDLDR, 18,26 

GPFAFT, 18 

IORTH,61,65 

IWRKSC, 23 

LATTDR, 18, 26 

NOTBIN, 18 

ORIBIN, 23 

PBCGST, 24 



PETRKP, 24 
SHPARi,SHPAR 2 ,..., 32 
TOL, 23 

ddscat.par - BETAMI,BETAMX,NBETA, 30 
ddscat.par - PHIMIN,PHIMAX,NPHI, 30 
ddscat.par - THETMI,THETMX,NTHETA, 30 
ddscat.par parameter file, 17 
diel.tab - see ddscat.par, 18 
Dielectric Frame (DF), 34, 42, 67 
dielectric function of target material, 26 
dielectric medium, 1 1 
dielectric tensor, 67 

differential scattering cross section dC sca /dfl, 64 
DOTORQ, 27 

effective radius a e g, 7, 19 

Electric field in or near the target, 68 

ELLIPSOID, 38 

error tolerance, 23 

ESELF, 71 

ETASCA, 27 

extinction efficiency factor Q cx t, 1 1 

FFT algorithm, 24 
FFTW, 24 
Fortran compiler 

optimization, 15 

GAMMA - see ddscat.par, 19 
GKDLDR, 26 

GKDLDR - see ddscat.par, 18 
GPFA, 24 
GPFAFT, 49 

GPFAFT - see ddscat.par, 18 
GPFAPACK, 71 

IDVERR, 17 

IDVOUT, 17 

infinite cylinder, 45 

infinite hexagonal column, 49 

IORTH, 65 

IORTH parameter, 6 1 

iterative algorithm, 23 

IWRKSC - see ddscat.par, 23 

IWRPOL, 68 

Lab Frame (LF), 28 
LATTDR, 26 

LATTDR - see ddscat.par, 18 
LFRAME, 55 

Magnetic field in or near the target, 68 
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memory requirements, 28 
Microsoft Windows, 13 
MKL 

Intel® MKL 
DFTI, 16 

MPI - Message Passing Interface, 16 

code execution, 2 1 
mtable - output file, 22 

Mueller matrix for infinite targets, periodic in 1-d, 66 
Mueller matrix for infinite targets, periodic in 2-d, 66 
Mueller matrix for scattering, 61, 63, 64 
MXNX,MXNY,MXNZ, 28 

NCOMP 

ddscat.par, 18 
NOTBIN - see ddscat.par, 18 

orientational averaging, 3 1 

nonrandomly-oriented targets, 3 1 

randomly-oriented targets, 3 1 
orientational sampling in j3, 6, and 30 

PBC = periodic boundary conditions, 57 
PBCGST algorithm, 24 

periodic boundary conditions, 44, Ab-A9, 51-53 

PETRKP algorithm, 24 

phase lag efficiency factor Q p h a , 1 1 

PIM package, 24 

polarizabilities 

GKDLDR, 26 

LATTDR, 26 
polarization - elliptical, 62 
polarization of incident radiation, 58 
polarization of scattered radiation, 64 
precision: single vs. double, 15 
PYD, 45, 48, 56, 57 
PZD, 45, 48, 57 

qtable - output file, 22 
qtable2 - output file, 22 
qtable2 file, 22 

radiative force efficiency vector Q pr , 27 
radiative torque efficiency vector Qr, 27 
READPOL, 68 

refractive index of target material, 26 
relative dielectric function, 1 1 
relative refractive index, 1 1 

scattering - angular averages, 58 
scattering angles 6 S , <fi s , 59 
scattering by tilted cube, 58 
scattering directions, 55 
scattering efficiency factor Q sca , 1 1 
SCATTERLIB, 71 



size parameter x = ka e g, 1 
source code (downloading of), 13 
Stokes vector (I, Q, U, V), 63 

Target Frame, 67 
Target Frame (TF), 28 
target generation, 32 

target generation: infinite periodic targets, 43 

target orientation, 28 

target routines: modifying, 43 

target shape options 

ANI_ELL_2, 35 

ANI_ELL_3, 35 

ANIELLIPS, 35 

ANIFILPBC, 44 

ANIFRMFIL, 34 

ANIRCTNGL, 35 

BISLINPBC, 44 

CONELLIPS, 36 

CYLINDER1, 36 

CYLNDRCAP, 37 

CYLNDRPBC, 44 

DSKBLYPBC, 46 

DSKRCTNGL, 37 

DSKRCTPBC, 47 

DW1996TAR, 38 

ELLIPSO_2, 38 

ELLIPSO.3, 38 

ELLIPSOID, 38 

FRMFILPBC, 43 

FROMJTLE, 32 

HEX_PRISM, 38 

HEXGONPBC, 48 

LAYRDSLAB, 39 

LYRSLBPBC, 49 

MLTBLOCKS, 39 

RCTGL_PBC, 51 

RCTGLBLK3, 39 

RCTGLPRSM, 39 

RECRECPBC, 51 

SLAB_HOLE, 40 

SLBHOLPBC, 52 

SPH_ANI_N, 41 

SPHERES.N, 40 

SPHRN.PBC, 53 

SPHROID_2, 41 

TETRAHDRN, 42 

TRILYRPBC, 54 

TRNGLPRSM, 42 

UNIAXICYL, 42 
target.out - output file, 43 
TFRAME, 55 

Vista, 13 
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wOOOrOOOkOOO.poln - output file, 22 
wOOOrOOOkOOO.sca- output file, 22, 23 
wOOOrOOOori.avg - output file, 23 
wOOOrOOOori.avg file - output file, 22 
Windows 7, 13 



